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This paper presents an overview of some techniques and concepts coming from dynamical system 
theory and used for the analysis of dynamical neural networks models. In a first section, we 
describe the dynamics of the neuron, starting from the Hodgkin-Huxloy description, which is 
somehow the canonical description for the "biological neuron" . We discuss some models reducing 
the Hodgkin-Huxley model to a two dimensional dynamical system, keeping one of the main feature 
of the neuron: its excitability. We present then examples of phase diagram and bifurcation analysis 
for the Hodgin-Huxley equations. Finally, we end this section by a dynamical system analysis for 
the nervous flux propagation along the axon. We then consider neuron couplings, with a brief 
description of synapses, synaptic plasticiy and learning, in a second section. We also briefly discuss 
the delicate issue of causal action from one neuron to another when complex feedback effects and 
non linear dynamics are involved. The third section presents the limit of weak coupling and 
the use of normal forms technics to handle this situation. We consider then several examples of 
recurrent models with different type of synaptic interactions (symmetric, cooperative, random). 
We introduce various techniques coming from statistical physics and dynamical systems theory. 
A last section is devoted to a detailed example of recurrent model where we go in deep in the 
analysis of the dynamics and discuss the effect of learning on the neuron dynamics. We also 
present recent methods allowing the analysis of the non linear effects of the neural dynamics on 
signal propagation and causal action. An appendix, presenting the main notions of dynamical 
systems theory useful for the comprehension of the chapter, has been added for the convenience 
of the reader. 
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I. INTRODUCTION. 

The present chapter aims to give an outlook of the various dynamical systems notions and techniques that are 
used while modeling Neural Network dynamics. Actually, there are a lot of such models. One reason is that there 
arc several levels of description and abstraction in this context : from a biologically realistic modeling of a neuron to 
neurons with a binary state; from an isolated neuron to Neural Networks, composed by several functional parts, each 
of them constituted by many neurons, and interacting in complex fashion, etc .... Another reason is that the Neural 
Network community is wide : from biologists, ncurophysiologists, pharmacologists, to mathematician, theoretical 
physicists, including engineers, computer scientists, robot designers, etc... Clearly the motivations and questions 
are different. Models that are designed to tackle a given problem may have very different structure and properties. 
It follows from these remarks that any attempt to "give an outlook of the various dynamical systems notions and 
techniques used while modeling Neural Network dynamics" is necessarily partial, biased et includes arbitrary and 
subjective choices. For sure, this chapter is subject to these restrictions. 

With this idea in mind we made the choice to explore the world of Neural Networks according to a specific map, 
represented in the Figures 1,2, localizing the various models studied in this chapter in a 3 dimensional space. We 
started from the obvious remark that a Neural Network is roughly made of neurons and synapses. But there are 
different levels of complexity and accuracy in the description of neurons and synapses. For neurons we used a model 
categorization along two axis. The first axis is relative to the proximity to biology. In this hierarchy, the Hodgkin- 
Huxley model is at the first rank (Section II. A). The Hodgkin-Huxley equations are derived in the section II. A 
and some aspects of their dynamical properties are briefly described in the sections II. C (examples of bifurcations 
occurring in the Hodgkin-Huxley model when a control parameter such as the external current is applied) and 

II. D (propagation of a spike along the axon). Before this, one remarks that the Hodgkin-Huxley equations can be 
reduced to a two dimensional dynamical system, taking various forms according to the modeling, but retaining in 
particular one of the main feature of the neuron: the property of excitability. The two dimensional excitable dynamical 
system obtained by reducing the Hodgkin-Huxley equations are easy to understand and provide fairly pedagogical 
examples. The excitable systems come therefore next in our hierarchy (section II. B). They allow one to capture some 
important dynamical aspects in neuronal behavior, such as spike generation, refractory period, threshold, and they 
exhibit various dynamical regimes observed in the experiments. After presenting the general structure of models for 
excitable membranes (section II.B.l) we discuss several canonical examples in neuron modeling. The first example is 
the Fitzhugh-Nagumo model (Section II. B. 2), then we briefly present the Morris-Lecar model (Section II. B. 3), and 
Integrate and Fire models (Section II. B. 4). 

These sections essentially deal with "spiking" neurons, namely the activity of the neuron is manifested by emission 
of action potential or "spikes" according to various pattern (individuals spikes, periodic spiking, bursting, etc . . . ). On 
biological grounds, this is certainly a fundamental aspect in neuronal dynamics. However, another description of the 
neuron can be made in terms of "firing rates" . The firing rate is the frequency of the spikes occurring during a certain 
time window of length T (typically, T ~ 100ms). It plays certainly also an important role in a certain number of 
neurological processes. For example, it is known since a long time (6), (7) that the firing rate of stretch receptor neurons 
in the muscles is related to the force applied to the muscle. However, during recent years, experimental evidences 
have suggested that this concept may be too simplistic to describe brain activity. It neglects indeed important aspects 
such as the information possibly contained in the exact timing of the spikes (1; 25; 91; 123; 125; 140; 147). Also, 
the reaction times in behavioral experiments are often too short to allow long temporal averages (see for example the 
experiments by S. Thorpe (151) on the vision). 

Nevertheless, firing rate models play an important role in the Neural Network community since they have been 
often used to model the collective activity of a neural assembly (8), (9), (44), and also to perform recognition tasks 
(90). Henceforth, we have included them in our table, and we have placed them after the spiking neurons in our rough 
hierarchy. In the examples described in the sections V,VI, corresponding to recurrent neural networks, the neuron 
is basically considered as an entity having an input and an output with a non linear transfer function (typically a 
sigmoid). This nonlinearity has several deep effects on the dynamics and a detailed example is described in section 
VI. 

Finally, if one makes the further approximation that the slope of the sigmoid function is infinite, one ends up 
with a binary state neuron (or Mac CuUogh and Pitts neuron (109)). Neural Networks with such binary "spin" like 
neurons had a great success (90) but we shall not discuss them in this chapter. 

The second axis of the table 1 takes into accoimt the collective aspects of Neural Networks. We establish a hierarchy 
ordering the models by increasing complexity in the neural population: one neuron, then a few neurons, then one 
population of weakly coupled neurons, then one population with arbitrary couplings (one could also consider several 
populations interacting with each others, but we do not consider this case in this chapter). 
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FIG. 1 Different levels of description of the neuron/network and techniques used to handle the dynamics of the models presented 
in this chapter. 



If one observes this space and asks which analytical methods allow us to describe the dynamics, one obtains the 
Table 1. A simple glance reveals that the methods discussed in this chapter essentially belong to three different 
domains of mathematics and physics: Dynamical systems theory, statistical mechanics and probability theory, and, 
at the intersection, ergodic theory. Also, one can remark that we essentially deal with the diagonal of this Table. 
As a matter of fact, when one moves away from the diagonal one meets, on one side, more and more trivial models 
(e.g. an isolated binary neuron), and, on the other side, more and more complex cases (a big population of many 
Hodgkin-Huxley neurons). In the first case, there is almost nothing non trivial to say, and in the second one, very 
little is known at least from the analytical point of view. In this chapter, we therefore choose some examples on the 
diagonal and we analyze the corresponding dynamics. 

There is actually, behind this choice, a fundamental aspect in modeling and analyzing Neural Networks, and 
more generally, modeling and analyzing the so-called "complex systems" . Complex systems are often composed by 
elementary units (in our case, neurons), having their own intrinsic characteristic dynamics and interacting with each 
others in a complicated way (nonlinear, non symmetric, with delays, etc . . . ). The intrinsic dynamics of the units 
can already be quite a bit complex (see, for example, the section II. C) so one may expect the collective dynamics 
to be even more complex. This is certainly true, but coupling the units give usually rise to a collective emergent 
behavior that one may characterize by the sentence: "The system as a whole is not reducible to the superposition 
of its elementary components" . This is usually due to non linear effects but this can also result from large numbers 
effects. Nevertheless, when one builds a dynamical system by coupling entities, each of them described by a lower 
dimensional dynamical system, the wisdom acquired when observing individual units is usually not sufficient to handle 
the collective behavior. The coupled system inherits characteristics that cannot be inferred from the knowledge of 
the uncoupled one. Also, some characteristics of the individual units may be hidden or may become irrelevant in the 
collective dynamics. These emergent effects can arise even if the coupling is weak. Starting from isolated neurons and 
"switching on" an interaction (synapses) between them, with an increasing intensity controlled by some parameter, 
the coupled system may, in some situations, exhibit a sharp, drastic change in its dynamics even if the parameter 
is small. This change usually corresponds to a bifurcation and it has often some analogies with phase transitions in 
statistical mechanics. Some prominent examples are presented in section IV. 
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The existence of emergence has two consequences. Firstly, this justifies somehow the simplifications inherent to 
modeling. If one desires to understand some emergent properties resulting from coupling neurons it might not be 
necessary to integrate all the features of the isolated neuron. It is often possible to drop some feature (preventing, for 
example, an analytic computation) and to capture nevertheless some important collective aspect. This outlines one 
important feature of the diagonal in table 1. When going from one "level of complexity" (detailed description of the 
neuron dynamics) to another level (coupling neurons) one often simplifies the characteristics of the neuron in order to 
have a tractable model. This is in some sense what we do when going from spiking Hodgkin-Huxlcy neurons to firing 
rate neurons. However another consequence results from the modeling process aiming to capture some characteristics 
considered as "relevant" and eliminate others considered as "details" . The mathematical structure and properties of 
the coupled model might be drastically different from the uncoupled one. This means that the tools, techniques or 
even philosophy adopted to handle the dynamics may change from one level of complexity to another. As we shall 
see, for example, the normal forms theory is quite a bit useful to handle dynamical changes in isolated neurons or in 
weakly coupled neural networks (provided some necessary assumptions are made) , but it is of little help in randomly 
coupled recurrent neural networks, at least before any prior treatment (such as the dynamic mean field equations of 
section VI. D). 

It results from these remarks that there is, currently, no general strategy to study Neural Networks dynamics. 
Nevertheless, as we shall see, dynamical systems theory, probability theory, statistical physics and ergodic theory can 
sometimes be used and combined to give partial solutions and can be tailored to build new tools and methods. A 
few examples are given in this chapter. 

Now, a few words about Table 2 below. It defines a third dimension in our classification space, where we define 
several levels of description for the "synapses" (interactions between neurons). The detailed physiology of the synapse 
is complex and, actually, there exists different types of synapses: chemical or electrical (gap junction). However, in 
most models the mathematical description is rough and, quite often, synapses are basically modeled in a way allowing 
to store information in the network, this information being extracted from the dynamical evolution of the neurons. 
Depending on the modeling chosen for the synapses, the dynamics can be very different, and their modification can 
induce drastic dynamical changes. In this chapter, we essentially give one example of the changes induced when one 
considers the different types of synapses presented in table 2, for recurrent networks (section V). We discuss first 
the convergence properties of the Cohen-Grossberg model when the synapses are symmetric (section V.B). Then 
we discuss the case of cooperative networks. The main result is a convergence theorem from Hirsch (87) which had 
recently some extensions in the field of genetic networks (74; 142). We also discuss in this section the notion of 
frustration resulting from the competition of excitatory/inhibitory effects. The section VI is devoted to the complete 
analysis of a recurrent model with asymmetric interactions, exhibiting complex regimes such as chaos. One can indeed 
go quite a bit deep in the description of the dynamics, by combining dynamical systems theory, statistical mechanics 
and ergodic theory (sections VI. A, VLB, VI. C, VI. D). This model exhibits interesting properties when submitted to 
Hebbian learning (section VI. E). We also present new developments characterizing the ability of such a network to 
transmit a signal. The basic tools is a linear response theory recently developed by Ruelle (131) (section VI. F). 
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FIG. 2 Different type of (formal) synapses considered in this chapter (section V). 



To conclude this introduction we would like to point out an important aspect. Many techniques described here 
have been developed out of the field of Neural Networks. But, in many cases they have been tailored or adapted to 
tackle specific problems in this field, and new methods have emerged. The interesting remark is that some of these 
techniques have now applications in other fields such as genetic networks, communication networks, or more generally 
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non linear dynamical systems on non regular graphs^, with a large number of degree of freedom (but finite). Some 
examples of applications to other fields are discussed in this chapter. 



II. SPIKING NEURONS AND EXCITABLE SYSTEMS. 



The activity of a neuron occurs by the emission of action potentials (or spikes) (see Fig. 3). In the simplest cases, 
they are controlled by ions (mainly Sodium (Na+) and Potassium (K-|-)) and their concentration around the nerve 
cell (see section II. A). An external stimulus causes Na-selective ion channel to open causing an influx of Sodium in 
the nerve cell. If the corresponding potential exceeds a threshold value (depolarization threshold) an action potential 
is generated. The action potential propagates then along the axon (section II. D). After the cell depolarizes, it must 
repolarize to its resting potential before it can depolarize again. This repolarization phase is controlled by an efflux of 
Potassium (repolarization phase) . This phase is followed by a refractory period where the neuron cannot be excited. 
The initial balance between Sodium and Potassium is restored by ionic pumps. Different models accounting for action 




FIG. 3 Typical action potential of a neuron. 



potential generation exist and some of them are described below. But, the core of all these models is certainly the 
Hodgkin-Huxley's that we describe in the next section. 



A. Hodgkin-Huxley neurons. 

The classical description of neuronal spiking dates back to Hodgkin and Huxley (89) . After extensive experimental 
studies these authors were able to propose a model for the dynamics of the giant axon of the squid. This constituted 
a signiflcant breakthrough in the description of action potential. At the time of their experiments (1952), the modern 
concept of ion-selective channels controlling the flow of current through the membrane was only one hypothesis 
among several competing others. Their model ruled out alternative ideas and gave correct predictive results of 
experiments that were not used in formulating the model. It reproduces and explain a remarkable range of data 
from squid axon, including the shape and propagation of the action potential, its sharp threshold, refractory period, 
anode-break excitation, accommodation and sub-threshold oscillations. Hodgkin & Huxley also proposed a set of 
equations modeling spikes propagation along the axon (see section II. D). They were in particular able to predict the 
propagation rate of spikes with a remarkable accuracy. The Hodgkin-Huxley modeling is generic, tractable and gave 
rise to new techniques and concepts. Consequently, the actual models of neural excitability are greatly influenced 
this work which resulted in a Nobel price (1961) for the authors. There is a large number of papers and books dealing 
with Hodgkin-Huxley model. Our main references are (48; 72; 85; 101; 104; 120) 

In their work, Hodgkin and Huxley start from the idea that the action potential results from transmembrane 
currents mainly constituted by Sodium (Na^) and Potassium (i^^) ions. Consider a neuron at rest in its natural 



^ In this way, the wisdom coming from the field of Neural Network is different (and complementary) from the knowledge acquired in 
parallel fields, such as coupled map lattices. 
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environment, namely in the intra cellular fluid where the Sodium and Chloride concentration is similar to sea water. 
One observes that, at rest, the Na~^ concentration is about 10 times higher outside tlic^ neuron than inside, while 
the concentration is about 5 times higher inside than outside. Assuming that the system is locally at thermal 
equilibrium with a temperature T, the difference in concentration between the inside and the outside, for the ionic 

species X, results in a potential difference Ex Vin[X] — Vout[X] called the Nernst potential and given by: 

where R = Afk = 8.315 J/ is the ideal gas constant, {Af = 6.02 x 10^^ is the Avogadro number, k = 1.38 x IQ-'^^J/K 
the Boltzmann constant), F = Me = 96500 C is the Faraday number (e = 1.602 x 10~'^^C is the charge of the 
proton), and [-''^lo^t (resp. [-'^Jj^) is the concentration of X outside (resp. inside) the neuron. With this convention, 
for positive ions, the effective electric force has the same direction as the force induced by the concentration gradient. 
For the giant axon of the skid and for a temperature T = 6.3 °C, the Nernst potential for Sodium and Potassium are 
respectively E^a ~ b6mV, Ek ~ — 77mV^. Moreover, taking into account the respective concentration of all ionic 
species the membrane potential is about —70mV at rest. Were the membrane to be permeable to ions, would one 
observe ionic currents through the membrane. These currents are not observed at rest, but arise during an action 
potential. Consequently, the ionic permeability of the membrane (conductance) depends on the neuron state (i.e. its 
membrane potential). 



In Hodgkin-Huxley modeling the (macroscopic) membrane conductances are determined by the combined effects 

of a large number of microscopic ionic channels located in the membrane. One considers a channel as an ensemble 
of independent gates (that can be of different type) with a binary open-closed state. Denote by pi = Pi{V) the 
probability that a a gate of type i is open. Then the conductivity Gx for channels of ionic species X, with gates of 
type i = 1 . . . iV, is proportional to the product of the probabilities pi that the gate i is open : Gx = gx YliLi Pi- 
where gx is the maximal conductance for channels of type X. Each pi depends on the potential V and on the fraction 
of open (pi) and closed (1 — Pi) gates. In the Hodgkin-Huxley model the time dependence of the Pi's is given by a 
master equation: 



— = a,{V)il - p,) - mp, = —^^^ (2) 

where at (resp. /Si) are the transition rates from close to open (resp. open to close) or gate inactivation (resp. 
activation). They have been empirically determined by Hodgkin and Huxley for each ion species. They are function 
of the membrane potential V (see eq. (13) below). In the second equality one introduces the natural quantities: 



ai{V) + /Ji{V) ai(V) + fii{V) 

where Ti is a characteristic time constant and is called the steady state activation. This is the value reached by pi 
when it is held at a potential V for a long period (say larger than the characteristic time r,). The solution of (2) is 
obviously: 



Piit) = PTiV) - iPTiV) - Pl)e-^ (4) 

Consequently, for a fixed V, pi has a simple exponential time dependence governed by r,. 

From their experiments Hodgkin-Huxley proposed to model the K conductance with an equation of the form: 



Gk = QKn^ (5) 

where gK is the maximum Potassium conductance. This corresponds to have a K channel with four independent 
gates of type n. The probability n is called the K activation variable. 
A similar equation can be written for the sodium: 



GjVo = QNam^h 



(6) 
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This corresponds to model a iVa+ channel with three gates of type "m" and one gate of type "/i" . m is the Na 
activation variable, and h is called the Na inactivation variable. The Na'^ ions can penetrate in the cell only if the 
m and h gates are both open (see Fig. 5). 



The membrane potential V is now given by Kirchhoff law 

iNa + Ik 



dV 



II 



(7) 



where I Ma, Ik are the sodium and potassium ionic currents through the cell membrane, II the leakage current (mainly 
composed by CI" ions) and I^xt is some external current (for example applied during an experiment). Cm is the 
membrane capacity {~ IfiF/cm^). The currents are given by the Ohm's law /, = Gi{V — Ei) where Ei is the Nernst 
potential of the species i = Na, K, L. 
Finally, the ionic currents are given by : 



C, 



dV 



dt 

1 dn 

1 dm 
W)~dt 
1 dh 



= -gNarrv'hiV - E^a) - 9Kn\V - Ek) - g^V - El) + 7ext 



a„(y)(l - n) - 0n{V)n 



amiV){l - m) - l3m{V)m = 
ahiV)il - h) - MV)h = 



TniV) 

m°°{V)-m 

Tm{V) 

h°°{V)-h 



(8) 
(9) 
(10) 

(11) 



The dynamical system (8-11) constitutes the complete Hodgkin- Huxley system. It involves a temperature dependent 
factor : 



(r-6.3) 



(12) 



This factor has the only effect of modifying the time constants in the equations for the activation/inactivation 
variables^. In the sequel we shall forget it and assume that the temperature is T = 6.3 °C ('y(T) ~ 1). 

The V dependence of the parameters an, f3n, am, Pm, cih, /3h was determined empirically by Hodgkin and Huxley. 
They found ^: 



am{V) = * 
ah{V) 



-(T/ + 45) 
10 

-(1^ + 60) 
10 



-(v+70) 

0.07e — 53 — ; 



-(V + 70) 

Pm{V) = Ae—T^ 



Pn{V) = 0.125e- 



-(V+70) 



1 



1 + e 



(-(V+40)) 



with: 



1 



if a; = 



(13) 

(14) 
(15) 

(16) 



In Fig. 4a have we drawn the time constants Tn,Th,Tm deduced from eq. (13) as functions of V, while in fig. 4 b 
the steady state values n°° , m°° , h°° as functions of V are shown. One notes in particular that the time constant for 
the Na activation variable is about one order of magnitude less than for the Na inactivation and the K activation, 
through the entire range. This means that the response in the m variable is quite a bit faster than the other variables. 
Consequently, during an action potential, when the voltage is high and m is large, it will take a while for h to decrease 



^ For a recent numerical work on the effects of temperature on the dynamics of a network composed by Hodgkin-Huxlcy neurons, coupled 
with gap junctions, see (155). 

^ In the literature one may find different forms for these equations depending on the zero of the potential. Here we have chosen it such 
that the membrane potential at rest is Vrest = — TOmF. One can also choose it such that y = at rest. 
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FIG. 4 Time constants T„,Th,r,n and steady state values n°°,m°°,h°° as functions of V. 



and for n to increase and contribute to the opposite K current. The mechanism of action potential emission is then 
the following. In the resting phase (a) the m, n gates are closed while the h gate is open. Therefore, sodium and 
potassium are neither leaving nor entering the cell (fig 5a). During depolarization, the m gates open fast allowing 
sodium to diffuse inside the cell, following the concentration gradient, while the n gates are still closed (fig 5b). This 
increases the membrane potential. Then n increases slowly, more and more K gates are open, generating an opposite 
K current. In the same time, h decreases and more and more h gates close, preventing sodium from coming into the 
cell (fig 5c). This corresponds to the repolarization phase. In the refractory period the m gates close, the h gates 
stay closed and the n gates stay open. It is not possible to excite the neuron in this phase (fig 5d). Finally, the h 
gates open, the ionic balance is restored by ionic pumps, and the resting state is once again achieved. If one draws 
the membrane potential versus time one obtains a picture similar to figure 3. The action potential is then propagated 
along the axon. The propagation equations are studied in section II. D. 




IN 

K+' 



FIG. 5 The various phase of the action potential in terms of the Hodgkin-Huxley equations. 

The preceding analysis is only qualitative but deeper mathematical investigations can be done (see section II. C) 
and numerical simulations can be performed. One observes spike generations but also periodic spiking, bursting etc... 
The Hodgkin-Huxley equations describe therefore the neural dynamics with a fantastic accuracy accounting of the 
wide variability in neuron activity. In particular, one predicts various situations observed in experiments. On the 
other hand they equations can be simplified giving rise to many models of formal neural networks. Despite this 
simplification (that can be quite a bit rough) it is still possible to obtain a huge quantity of information about the 
neural dynamics. In the next section we present a few models derived from Hodgkin-Huxley equation and capturing 
one of the main feature of the biological neuron: excitability. 
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B. Reducing the Hodgkin-Huxley equations. 

1. General structure of excitable membrane. 

Most models for excitable membrane retain the general Hodgkin-Huxley structure (eq. (8)-(ll)) and can be written 
in the form. 



dV 
dt 



C'm-j^ — —IjoniYi X\, ■ . . , Xn) + lext — ~ IkjYi ^1; ■ • ■ ; ^w) + ^ext^ (17) 

fe=l 

h = 9kC7k{V,Xu...,Xr,){V -Ek), k = l...N, (18) 



dt Ti{V) 



(19) 



where V denotes membrane potential, the membrane capacity, li^n is the sum of ionic currents, I^xt external 
or applied current. The variables pi are used to describe the fraction of open channels of type i. is the characteristic 
time that the ions of type i need to reach the rest state pfiV). In the Ohm's law (18), Ik is the current for the k 
th ion species, gk is the maximal conductivity for the ions channels of type k, Gk is the product of gate fc-channels 
activity, and Ek is the Nernst equilibrium potential. In some situations it is fundamental to have an accurate models 
of the neuron excitability, if one seeks, for example, to account for rather detailed aspects of spike shape, dependence 
upon many pharmacological agents, etc .... However, in many cases a rough description is enough to capture the main 
qualitative and quantitative aspects of the dynamics of excitability. Consequently, one can reduce the complexity 
of the set of equations (17,18, 19) in order to obtain an analytically tractable model. Henceforth, many models of 
neuronal dynamics are reduction of these general equations. 



2. The FitzHugh-Nagumo model 

In this spirit FitzHugh (65) and independently Nagumo, Arimoto et Yoshizawa (119), considered reductions of the 
Hodgkin-Huxley model and introduced an analytically tractable two variables model. 

The basic observation is the time scale separation between the variables V,m,n,h in eq. ((8)-(ll)). According to 
Fig. 4 the characteristic time for Sodium activation is so fast compared to the other variables that one may consider 
m essentially as a constant. This eliminates the variable m. Also, FitzHugh observed that h + n \s essentially a 
constant ^ 0.8 during the action potential. Consequently, one can eliminate one more variable. One finally obtains a 
model of the form (for the detailed reduction see e.g. (4), (102), (101), (72), (126)): 

e- = hiv,w) (20) 
dw , , 

= 9x{v,w) (21) 

C 

where e = — r- is typically small. The index A refers to the control parameters of the system. In the FitzHugh- 

maxy Tn(V) 

Nagumo model f\{v, w) = v — v'^ — w + / is a cubic polynomial in v and is linear in w, while gx{v, w) = {v — a — bw). 
The parameters A = (a, b, I) are deduced from the physiological characteristics of the neuron. It can also be useful to 
consider the dynamical system 

^ = fxiv,w) (22) 
dw , , 

— = tgx{vM (23) 

obtained from (20,21) by a time rescaling t ^. 

The system of equations (20,21) is the canonical form for excitable systems. That is why we used the "generic" 
variables namely v, w instead of V, n. They are usually called excitation and recovery variables. The excitation 
variable governs the rise to the excited state while the recovery variable causes the return to the steady state. Since 
e is typically a small parameter, there is a separation of time scales between the two variables. 
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On technical grounds, the analysis is simplified by the two dimensional geometry of the phase space. Indeed, in 
the phase plane, the slope of the trajectory of a given point is ^ = ^^^^'^j and consequently the phase portrait can 
easily been drawn. In particular a trajectory is vertical (resp. horizontal) at the points such that f\{v,'w) = (resp. 

gx{v, w) =0). The set of points Ny =* {{v, w) \ fx{v, w) = 0} (resp. =^ {{v, w) \ g\{v, w) = 0}) is composed by a 
union of curves called the v-nullclines (resp. w-nullclines) . Thus, the fixed points of (20,21) are at the intersection 
of nuUclines. More generally, the shape of the nuUclincs gives strong informations on the dynamics. As shown below 
the nuUclines shape changes when the parameters A are varying, leading to bifurcations for some values of A. 

When e is small one uses an additional property to analyze the dynamical system (20,21). Setting e = in (20,21) 
one obtains fx{v, w) = 0; ^ = g\{v, w). This means that, whenever it is possible, v is adjusted rapidly to maintain a 
pseudo-equilibrium corresponding to f{v,w) = and plays the role of an implicit parameter in the evolution of w. In 
other words, the point {v, w) moves slowly along the (stable) branches of the v nuUclines. These branches compose 
the so-called "slow manifold" : it is only "on" (or very near) this curve that the motion of the solution curves is not 
very fast in a nearly horizontal direction (sec e.g. Fig. 6). 

On the other hand, away from the Nw nuUcline, the vector field is essentially horizontal and one has a fast motion 
of V. Indeed, a time rescaling t ^ gives the system (22,23). Then, setting e = one can approximate the (regular) 
trajectories of the system (20,21) by the (non regular) trajectories of the degenerated system: 

^ = fx{vM (24) 
dw „ 

^ = ' ^''^ 

where the vector field is horizontal with a norm f\{v, w). 

The trajectories of the real system arc composed by pieces coming from these two approximations. There arc 
theorems controlling how the real trajectories of (20,21) are close to the piecewise trajectories, for a sufficiently small 
e allowing to obtain the characteristic trajectories of the initial system from the solutions of the degenerated system. 
This is the essence of the singular perturbation theory developed by Mischenko & Rozov (118). 

To illustrate this, let us start we a simple example used as a preliminary step to analyze later on the FitzHugh- 
Nagumo equations: 

e— = v-v^-w (26) 



dt 
dw 



v-a (27) 



The i;-nullcline is given hy w = v — while the w-nuUcline is the vertical line v = a. The nuUclines and the flow of 
(26) are depicted fig. 6. Due to the smallness of the parameter e, the flow is essentially horizontal^ (^ ^ 0) except 

close to the w-nuUcline. Crossing the f-nuUcline (resp. the w-nuUcline) makes the v component of the flow (resp. 
the w component) changing its sign. The v nuUcline has two "stable" branches denoted by N^. Namely the flow is 
attracted in a neighborhood of these branc;lics and stays a long time in this neighborhood, moving slowly upward for 
the + branch and downward for the branch — . In the case of the + branch the flow finally reaches the extremum. 
Then it moves fast to the other branch. The middle branch is called the unstable branch. As discussed below it acts 
(roughly) as a threshold for spike generation. 

The point A = (jja = a- Wj\ = — a + a"^), where the nuUclines intersect, is a fixed point. The eigenvalues of the 

corresponding Jacobian matrix DFa are Ai^2 = 3a =*=V(^ ) Consequently, the eigenvalues are complex 

for a €] ^ — "^^"^^ [U] ^j^*^"* , ^'^'^^'^^ [ and real otherwise. Moreover, A is stable when \a\ > and 

unstable otherwise. More precisely, this is a sink (Ai,A2 < 0) for o e] — oo, — U +oo[, a stable 

focus (5R(Ai,2) < 0) for a e] - ^i±|^, -^[U]^, ^J^^^[, a center (5R(Ai,2) = 0) for \a\ = ^, an unstable focus 



Note that, strictly speaking, the vector field on the y axis is not horizontal (its component are (— a)). However, for simplicity, we 
are drawn it horizontal, assuming a very small e value. 
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FIG. 6 NuUclines and vector field for the toy model (26). This a qualitative drawing and the phase portrait has been drawn "by 
hand". Consequently, the arrows representing the vector field are drawn as indicators. The picture is not scaled. In particular, 
the vicinity of the slow manifolds (in green) is of order e. Practically, the trajectories near the slow manifold can essentially be 
considered as being "on" the slow manifold. 



(5R(Ai,2) > 0) for a e] - -^i^[u]i/i^, ^[, and a source (Ai,A2 > 0) for a G h V^^' V^^l (^^^^ 
the appendix for more details about the classification of fixed points). 

Assume now that we are in the situation depicted in Fig. 7a, with a < — The system is at rest in A. Now, 
we excite it moving A to B = {vb,wb- There are two possibihties. Either wb > 7^, then the excitation relaxes 

3\/ (3) 

down to the rest state (Fig. 7a). Or wb < 7^. Then we have the situation depicted in fig. 7b. The trajectory 

3V (3) 

flows rapidly parallel to v until it approaches the u-nuUcline and crosses it in C. Then it follows slowly the stable 
branch {C,D). At this point, the v flow is zero while the w flow is positive. Consequently, the trajectory leaves the 
nuUclincs, and is fast driven by the flow until the point E. It follows then the stable branch {E, A) until the rest state 
A. The corresponding trajectory of v is depicted in the inset of Fig. 7b. It has a spike shape where one recognizes 
the equivalent of the depolarizing phase {B,C), the repolarizing phase (C, E), and the refractory period {E, A) of the 
flgure 3. Consequently, this simplified model gives already a fairly good example of an excitable dynamical system. 

Note that the dynamical system (the neuron) is more sensitive to excitation when the fixed point A is closer to the 
local extremum Mi = (— ~ 3^^7(3) nuUcline (resp. M2 = (^, ^^^^ )), namely when the control parameter 

a is close to the bifurcation value a = (resp. a =^ "^)- ^^^^ way, one may consider that excitable neurons are 
dynamical systems close to a bifurcation point. This idea is further developed in section IV. This dynamical system 
has moreover an additional feature which makes it relevant to neuronal dynamics. Assume now that \a\ < Then 
the rest state A is unstable. If we slightly perturb A one generates a periodic activity depicted in fig. 7c. 

For general systems of the form (20,21) the nuUclines have a more complex shape and the dynamics is richer. It is an 
interesting exercise, illustrating the spirit of dynamical systems theory, to start from the system (26), and to ask what 
are the changes induced in the dynamics by deformations of the nuUclines. Let us do this for the FitzHugh-Nagumo 
model. 

^ = v-v^'-w + I (28) 
at 

dw , , , 

— = e(v-a-bw) (29 
at 

It is deduced from the system (26) by translating the u-nuUclines with a vertical displacement / and by tilting the w 
nuUclines which becomes the straight line w = for 6^0. From a qualitative point of view one can figure out 
without any computation which type of novelties will be induced by these changes. As shown in Fig. 8, 9 we can for 
example have appearance/coalescence of pairs of fixed points by saddle-node bifurcations and bistability. 
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FIG. 7 Examples of possible behaviors for the equation (26) in response to a perturbation of the rest state A. Fig. 7a. 
Relaxation to the rest state A. Fig. 7b. Spike emission. Fig. 7c. Periodic spikes train emission. 



4. 



Nw 



BC„* 



..Nw 



FIG. 8 Saddle node bifurcation and bistability in the FitzHugh-Nagumo model (28) when the parameter b increases. Note that 
the slope of the w nuUcline is j-. The same remarks as in Fig. 6 holds for the scaling of the arrows. 



On more general biological grounds, and though the FitzHugh-Nagumo equations are a simplification of the 
Hodgkin-Huxley equations, they exhibit some typical behavior of the real neuron. Let us list a few examples. 

• Action potential emission and threshold. The first observation is that a suitable input current can generate 
an action potential. Consider the case depicted in Fig. 10. There is a unique stable fixed point A. Consider 
now the line labeled by S. This line is called the threshold separatrix since it separates solution curves 
that represent action potentials from those that do not represent action potentials (48). This curve is not 
sharply defined here (see the discussion of type I excitability for a definition) but it is very close to the 
unstable branch and, between the minimum and the maximum of the v nuUclines, it essentially corresponds 
to the set {{x,y) \ fv{v,'w) = fw{v,w)}, where the vector field makes an angle of 45 ° with the v axis. Let 
us now consider the situations corresponding to the case 1 and 2 in Fig. 10. One perturbs the rest state 
by changing the membrane potential such that v is close to 5, but in the case 1 the perturbed point is 
"above" S and in the case 2 it is "below" S. Even if these two points are close to each other, the vector 
fields have a different orientation since the angle of the vector field with the v axis is, in the case 1 larger 
than 45 ° and in the case 2 it is smaller. This has the following consequence. In the case 1 the neuron 
returns to equilibrium without emitting a spike. On the other hand, in the case 2 the trajectory has to 
make a big excursion before returning to the rest state: there is a spike emission. The horizontal distance 
from A to S corresponds therefore to a threshold value 6. Note however that the concept of threshold, 
corresponding to a sharp transition, is questionable, in the Hodgkin-Huxley model, since there is no real clear 
cut firing threshold (see (105; 127); see also the discussion below about type I and type II models of excitability). 



• Existence of a refractory period. The Figure 10 also exhibits two regions labeled by AR for "Absolute 
Refractory" , and RR for "Relative Refractory" . These regions are defined as follows. Assume that the neuron 
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FIG. 9 Bifurcation diagram corresponding to Fig. 8. X corresponds to the projection on the A^'j, nullcline. bo is the critical 
point. 




t 



FIG. 10 Spike emission in the FitzHugh-Nagumo model. 



is spiking. If the corresponding point in the phase space is in the region AR, any further positive increase in 
the membrane potential will not be able to generate a new spike. On the other hand, in the region RR a a 
spike can be generated provided the clamped potential is strong enough. 

• Anodal break excitation. Assume that an action potential is generated and, during this, an external potential 
(anodal shock) is applied at the instant where the system is the point P in Fig. 11, with the effect to move P to 
P' . If the shock is large enough such that P' is on the left of the threshold separatrix, the action potential is abol- 
ished by the anodal shock. This phenomenon has been observed experimentally (see (48) and references therein). 

• Spike emission by hyperpolarization. Assume now that we apply a negative current / < in the situation where 
the system is initially at rest, with a stable fixed point A (Fig. 12). The cubic nullcline moves downward and 
A moves to A' . If we removes the current, the cubic moves upward. But then A' is no more a fixed point. Its 
trajectory is described in Fig. 12. This corresponds to a spike emission. 

• Periodic sequences of action potential. Assume now that we apply a positive current I > 0. For sufficiently high 
I A becomes unstable. Then the slightest excitation generates a periodic emission of spikes. It is indeed possible 
to show rigorously, using Mischenko and Rozov theorems combined with the Poincare-Bendixon theorem (78), 
that there exists a stable limit cycle (depicted Fig. 13). 

What happens now if we go on deforming the nuUclines ? For example, one can bend the line corresponding to the 
w nullcline transforming it into a parabola: this is the deformation of lowest non linear order. It is quite interesting 
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FIG. 11 Anodal break excitation in the FitzHugh-Nagumo model. 
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FIG. 12 Spike emission by hyperpolarization in the FitzHugh-Nagumo model. 



to remark that this leads to a system exhibiting neural excitability^ of type I and II. Indeed, the response of a neuron 
to permanent current stimulus can generate a periodic train of spikes with a determined frequency. In this case, one 
distinguishes two types of such excitability (this classification was proposed by Hodgkin in 1948). 

• Type I excitability. The spike train is generated with an arbitrary small frequency, depending on the applied 
current (Fig. 14). From a dynamical point of view, such type of excitability can be generated by the scenario 
depicted Fig. 15a, b,c. The variation of a control parameter (here the applied current) moves the v nuUcline 
such that a saddle-node bifurcation on a limit cycle occurs. For a critical value I — Ic there is an homoclinic 
connexion on the fixed point A. Consequently, the period is infinite (and the frequency is zero). Note that the 
amplitude of the cycle is independent of /. 

In figure 15a we have also qualitatively plotted the separatrix S which is here the stable manifold of B. Clearly, 
a perturbation to the left of S does not generate a spike, while a perturbation to the right corresponds to a 
trajectory making a big excursion around the unstable fixed point C, before returning to the rest state: this 
corresponds to a spike. 



^ Note that type II excitabihty exists already in the previous case. 
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FIG. 13 Periodic sequences of spikes in the FitzHugh-Nagumo model. 



FIG. 14 Variation of the spikes frequency with the control parameter (applied current in Fig. 15a,b,c, 16a,b,c). Fig. 14a. Type 
I excitability. ) Fig. 14b. Type II excitability.). 



• Type II excitability. The spike train is generated with a frequency staying a specific domain (Fig. 14b). From 
a dynamical point of view, such type of excitability can be generated by the scenario depicted Fig. 16. The 
variation of the applied current moves the x nuUcline such that a Hopf bifurcation occurs. The frequency 
depends slightly on / and the amplitude increases like the square root of the parameter distance to the critical 
value, as long as one stays close to the bifurcation point. Note that the example depicted in Fig. 16 does not 
use the fact that Ny has a quadratic shape. Actually, the same is obtained with a straight line. 



3. The Morris Lecar model 



The previous examples may look quite abstract since we deformed the nuUclines freely, without paying much 
attention to the biological relevance of this operation. Actually, there exist biologically plausible models exhibiting 
the behaviors presented above. An example is the Morris and Lecar model (62; 117) which was formulated in the 
context of the electrical activity of the barnacle muscle fiber. The Sodium channels are replaced by Calcium channels. 
One calls m the activation variable. The Calcium conductance is given by Gca ~ gcan^iV). There is no inactivation 
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FIG. 15 Type I excitability. Schematic example of a dynamical system exhibiting type I excitability. 
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FIG. 16 Schematic example of a dynamical system exhibiting type II excitability. 



variable h. The dynamics is given by: 



where : 



C ^ 

dt 

dw 

m 



-gCamoo{V){V - Eca) - gKw{V - Ek) - .9l(^ -El)+I 



Woo{V) 



1 + tanh 

1 + tanh 
1 



V -Vi 
Va 



cosh ( 



(30) 
(31) 

(32) 

(33) 
(34) 

(35) 



w is the fraction of open channels. This set of equations as a large number of parameter that one may vary in 
order to study the behavior of the neuron when physical characteristics, such has Vi, V27 are varying. However, 

from an experimentalist point of view, the only free parameter is the external current /. 

The V nuUcline corresponds to a situation where the applied current exactly cancels the ionic current. It is given 
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by 

/ = 9cawUy){V - Eca) - gKwiV - Ek) - g^V - El) 

It has a cubic shape and a variation of / as simply the effect of translating it paraUel to the V axis. The w nuUchne is 
the activation curve w = WooiY). This model displays a wide variety of dynamics such as spikes, oscillations emerging 
with zero or non-zero frequency and bistability. 



4. Integrate and Fire models. 

A convenient and simple model producing spikes is the so called leaky integrate and fire model. Consider the circuit 
drawn in Fig. 17. The device D is conducting when the potential is above a threshold 9 and has an infinite resistance 
otherwise. It acts therefore as a potential dependent switch. The total current is I ^ I n + Ic — + ^llt- Using the 
time constant = RC one obtains the equation of the leaky integrate and fire model: 

du 

Tm-j- = ~uit) + Rlit) (36) 
dt 

with the additional condition that u cannot increase above 9. Starting, say, from a zero potential u, u{t) increases 
until it reaches the threshold value 9. Then D switches on and the capacity unloads. Consequently, the potential u 
decreases exponentially fast, u is interpreted as a membrane potential and as the membrane time constant of the 
neuron. 
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FIG. 17 Schematic circuit of the integrate and fire model. 



In integrate and fire models, the form of the action potential is not explicitly described. Instead, one models the 
situation above by saying that, when the potential u reach the value 9, at some time tf, it is immediately reset to a 

new value Ur u{t'^) < 9 while a spike is emitted. Then, the membrane potential keeps the value Ur for a time 
corresponding to the refractory period. In this sense, spikes are formal events characterized by the "firing time" tf. 
A more general version of (36) is a non linear integrate and fire model (5): 

rm^^Fiuit))+G{uit))Iit) (37) 

where F, G are non linear functions of u. 

Though the integrate and fire model is a rough modeling of a spiking neuron it has several advantages. Firstly, the 
linear model (36) is exactly solvable. The potential u{t) resulting from an excitation with a time dependent current 
is easily found. For example, the current after a spike arising at time ti and before the next spike (u(t2) = 9) is given 
by: 

u{t) = Ure^^*-^^ + ^ [ e("^)/(t- s)ds; <e[ti,t2[ 
^ Jo 

Also, it is easy to model a network of integrate and fire neurons ^. In this framework the neuron i receives the 
spikes coming from other neurons, and the total current Ii(t) is the sum of spikes coming from each neuron j weighted 



Note however that the equation (38) holds in a more general setting. 
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by a quantity Jy roughly modeling the synaptic connexion between j and i: 

ii{t) = Y,j,, Yl «(^-*"0')) (38) 

where t„{j) is the n-th time of firing of the neuron j, a is a function modeling the spike, and the sum X]nQ)=i^ 
corresponds to an integration over a small time window. The spike function a can have different forms, but the 
simplest one is a Dirac 5 distribution, corresponding to have an "instantaneous" spike. 

Note that an equation with the form ( 38) is particularly well suited for a stochastic approach, where the firing times 
are randomly distributed e.g. according to a Poisson process. An example of this is given in Chapter II. Most of the 
analysis use a stochastic approach. However the evolution can also be investigated in a deterministic context, where 
the firing condition is determined e.g. by an Heaviside function. Then one has to handle a deterministic dynamical 
system with singularities. 



C. Qualitative analysis of the Hodgkin-Huxley equations. 

We now return back to the Hodgkin-Huxley equations. The analysis made in section II. A was only quantitative. But 
it has allowed us to understand the spike generation, by simple arguments on the characteristic times of the variables 
m, n, h, and their interpretation in terms of probabilities that a gate of a given ionic species is open or closed. A further 
analysis requires however to consider the complete non linear dynamical system (8-11) and its dependence in control 
parameters such as the external current /. Actually, the simplifications made in section II. B lead us to find several 
situations having a correspondence with experiments on real neurons. Since the equations (28) are a simplification 
of the Hodgkin-Huxley system, one expects to observe similar effects in the dynamical system (8-11). However, the 
reduced systems were two dimensional while the Hodgkin-Huxley system has four dimensions. Therefore, bifurcations 
and dynamical regimes (such as chaos) occurring in phase space having more than two dimensions are not observed 
in reduced systems like (28). Rinzel and Miller (128) first gave evidence of this. Doi and Kumagai (55) recently 
showed the existence of chaotic attractors in a modified Hodgkin-Huxley model that changes the time constant of one 
of the current by a factor 100, and, more recently, Guckenheimcr and Oliva (80) showed rigorously the existence of a 
Smale horseshoe (hence of chaos) in the Hodgkin-Huxley model with its original parameters. Finally, the reduction 
performed to obtain the equations (28) used several simplifications that can be discussed and that may bring some 
exogenous properties, not present in the initial model. 

For all these reasons, there is a clear need to perform an analysis of the Hodgkin-Huxley system. Obviously, it is 
always possible to make numerical simulations of this dynamical system and many papers have been written on the 
subject (sec for example (113) and reference therein). Also, there exists currently a lot of "on line" simulators on 
the Internet (71; 86). However, analytical results are also useful since they allow in particular to locate bifurcations 
points. This is certainly useful because this permits to reduce the explored area in the (huge) parameters space and 
to locate small regions that could be missed by a discrete sampling in a numerical simulation. In this section we 
present one example of such an investigation, due to Guckenheimer & Labouriau (79). This paper presents actually 
an approach combining rigorous methods from dynamical systems theory with numerical tools of formal calculus (for 
more details on this type of approach see also (81)). This allows the authors to draw a bifurcation diagram in a two 
dimensional parameter space corresponding to the potassium reversal potential^ vk = Vrest — Ek and to the current 
/ (the reversal potential of Sodium and Potassium can indeed be controlled experimentally (88; 97)). Consequently, 
the bifurcations presented are generic codimension one and two bifurcations. Actually, the bifurcation diagram 
presented in Fig. 18 presents an overwhelming richness of dynamical behaviors in a rather small parameter space 
region. This is a "zoo" in which one meets basically all species described in standard textbooks about bifurcation 
theory (78; 129) (see the appendix) plus some more "exotic" individuals such as the twisted saddle loop bifurcations. 
This is one reason why we have chosen this example: it shows how deep the dynamical systems analysis can go and 
how rich are the Hodgkin-Huxley equations. Additional references are (55; 82; 107). 

Let us start from elementary remarks. It is easy to show that the asymptotic solutions of eq. (8)-(9) are contained 
in the set {m, /i, nG [0,1]^ x i/ G — r, i/+ -t- r}, for some r > 0, and where = min(r/ivai i^K, J^l) and = 



In the paper, the variable v corresponds to a clamped potential with the opposite convention as in the section II. A (see note 3) 
V = Vrest — V where V is the membrane potential. 
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max(z/jva, j/x, j^l). Fortunately m,h,n stay dynamically in [0,1]"^ (these are probabilities !!). Indeed, if ra (resp. 
h,n) is equal to zero ^ > and if m = 1, ^ < 0. Also, \{ v > v+, ^ < Q and \i v < v^, ^ > 0. As 

t — > oo, 'm,,h,n — > mrx,,hao,noo- Consequently, if v* is a equilibrium of eq. (8) then [v* ^moo{i^*),h^{h'*),noo{v*) 

is an equilibrium of eq. ((8)-(ll)). Also, ^ = ^ G{v* ,moo{i^*),hoo{i^*),n^{i^*)) '= /(i^*) = I- Consequently, 
there exists a unique value of / for which A = {v* ,mao{v*),hoo{i^*),noa{i^*) is an equilibrium. When vk has the 
value found by Hodgkin-Huxley, / is monotonic and (8-9) has a unique equilibrium for each value of /. For fixed 
lower values of vk there are two saddle node bifurcations as / is varied, creating a region with three equilibria and 
corresponding to multistability (as in the example depicted in the previous section, Fig. 8). The two curves of saddle 
node terminate at a cusp point. These curves are obtained by varying v* , considered as a parameter and taking into 
account the transversality conditions TSN1,TSN2 in the appendix. In particular the determinant of the Jacobian 
matrix has to vanish. Given the equilibriiuii point, one also obtains the parameters value where Hopf bifurcations 
occur. Hopf bifurcation requires that two complex conjugate eigenvalues appear or disappear. This corresponds 
to conditions on the coefScient of the characteristic polynomial of the Jacobian matrix. This polynomial has the 
form + cs-T'"' + C2X^ + cix + cq. Considering v* as a control parameter and solving simultaneously the fixed point 
equations and the transversality conditions (TH1,TH2) in the appendix) one finds the set of parameters vk^I where 
a Hopf bifurcation occurs. At the intersection of the Hopf bifurcation line and the saddle-node bifurcation line, a 
Bogdanov-Takens bifurcation occurs (see the appendix). One observes also global bifurcations: collapse of two limit 
cycles, homoclinic connexion at the Bogdanov-Takens point, twisted saddle loop, degenerate Hopf bifurcation, etc ... 
The various bifurcations are depicted in Fig. 18. We used the following nomenclature (from (79)). For a description 
of the corresponding bifurcations see the appendix. 

Codimension one bifurcations. 

• sn: Saddle- node bifurcation: two fixed point coalesce and disappear (resp. appear), see Fig. 46 in the appendix. 

• h: Hopf bifurcation. A fixed point changes its stability and a limit cycle appear with a radius increasing with 
the control parameter (resp. a limit cycle decreases until it is reduced to a point and disappear while the point 
at the center changes its stability), see Fig 49 in the appendix. As discussed above this corresponds to type II 
excitability. 

• si: Saddle-loop or homoclinic bifurcation. The amplitude of a periodic orbit increases until it captures a saddle 
point and disappears, its period tending to infinity when the control parameter tends to the critical value. As 
discussed above this corresponds to type I excitability. 

• tsl: Twisted saddle-loop bifurcation. In dimension larger than two an orientation reversal along a homoclinic 
may occur. The homoclinic orbit is a two dimensional ribbon which is invariant under the flow with tangents 
in the directions of the weakest contraction at the saddle point. A twisted saddle loop occurs if the ribbon 
is not orientable.This bifurcation is also met in physical experiments about Rayleigh-Benard convection in a 
small geometry (see (94) for a mathematical analysis). Note that this bifurcation is usually related to period 
doubling. Also, for any n value, n integer, there exists a dynamical system, arbitrary close to the bifurcating 
system, having homoclinic connexions with loops of order n (see (95)). The dynamics can therefore be quite 
complex in the vicinity of this bifurcation. 

• dc: Double cycle or saddle-node bifurcation of cycles. Two periodic orbit coalesce and disappear. 

• pd: Period doubling bifurcation. A periodic orbit changes its stability, while a periodic orbit of twice its period 
coalesce with the bifurcating periodic orbit. 

Codimension two bifurcations. 

• c: Cusp. Three equilibria coalesce into one (see Fig. 46 in the appendix). 

• tb: Takens-Bogdanov bifurcation (see appendix. Fig. 50). 
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FIG. 18 Bifurcation diagram of the Hodgkin-Huxley equations when varying the parameters /, vk- This figure has been drawn 
"by hand" from the Figure 1 in (79). Stable equilibrium points are shown as black dots, unstable focus as white dots, stable limit 
cycles are closed curves with solid lines and unstable periodic orbits are dashed lines. One dimensional unstable manifolds of 
equilibrium points are shown together with curves of the "weak stable manifolds" of equilibrium points with three dimensional 
stable manifolds (see e.g. in the "tsl" and "pd" regions). 



• nsl: Neutral saddle-loop or homoclinic bifurcation. A periodic orbit changes its stability in a saddle loop at a 
point where the sum of the eigenvalues of the Jacobian matrix is zero. 

• tnsl: Twisted neutral saddle-loop bifurcation. 

• snl: Saddle-node loop. 

• dh: Degenerate Hopf bifurcation. 

We have also represented some qualitative changes in the dynamics arising when varying the Nernst potential Vk 
in Fig. 20 a,b,c. 

These results illustrate the complexity of the dynamics occurring in the Hodgkin-Huxley equations. There are many 
possibilities for the spiking patterns when the parameters are changed. One may however ask about the biological 
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FIG. 19 (a) Untwisted saddle loop, (b) Twisted saddle loop.). 




FIG. 20 Bifurcations occurring when following the paths I (Fig. 20 a), II (Fig. 20 b) drawn in figure 18. The corresponding values 
for the potential Vk range from —5.155 to —5.129 mV. 



relevance of these results. Note that in Fig. 18 the usual value of Vk ~ lOmV is far on the right of the graph and 
does not appear in a scaled figure. Indeed, the region corresponding to the path II ranges from —5.155 to —5.129 
mV. Thus its width is of order 20 fj,V ... and the potential is negative... Thus, some of these regimes may be difficult 
to find experimentally, since they correspond to very tiny regions in the parameters space and quite unusual value of 
parameters^. Another related question is: what happens when coupling such neurons ? For example, do the regions 
sl,pd,tsl, exhibiting a complex behavior, still exist when considering a neural network of Hodgkin-Huxley neurons ? 
We shall see in this chapter that coupling neurons with complex dynamics does not necessarily imply that the coupled 
dynamical system will have a complex dynamics. On the opposite, coupling neuron models with a simple evolution 
may lead to a complex evolution. 



Note that the Authors of (79) also explored the changes induced by a variation of the Potassium conductance gj^ but we do not discuss 
this here. 
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D. Axon propagation. 

The Hodgkin-Huxley equations (8-11) describe the behavior of a smah piece of neuron membrane. From the fun- 
damental laws of Physics, one can use them to obtain an equation describing the propagation of the action potential 
along the axon. One can in particular obtain the propagation speed. In this section we derive the propagation 
equation. We then discuss the existence of propagating solutions in a simplified version of the propagation equations, 
based on the FitzHugh-Nagumo model. 

Let V be the local membrane potential and R the resistance per unit length (as discussed in section II. A it 
depends on V). For simplicity, we shall use in this section the convention where F = at rest and we shall set 
Vx = Vrest — Ex where X = Na, K, L and Ex is the Nernst potential. Denote by x the coordinate longitudinal 
to the axon. One decomposes the current in the membrane into an longitudinal part (ia) and a transverse part i^- 
From local charge conservation one has: ia{x + dx) = ia{x) — im{x) ^ = —imix), while the Ohm's law writes: 
V{x + dx) — V{x) = —Ria{x) => ^ = —Ria{x). Consequently: 

= Rim{x) (39) 
The local transmembrane current is given by the Hodgkin-Huxley system (8-11): 



dV dV 

imdx = dx{Cm-^ + lion) = ^rn^dx + S{x) [gNam^h{V - Vno) + QKu'^iV - Vk) + gdV - Vl)] (40) 



where S{x) = 2TTr{x)dx is the membrane surface per unit length and r{x) the axon radius at x. Finally, the equations 
describing the spike propagation along the axon are: 

" Cm^ + 27:rix)[gNam^hiV-VNa) + 9Kn\V-VK)+9L{V-VL)] (41) 
- =My)(l-n)-/3„(y)n=-lJ^ (42) 

^ = ..(y)(i - m) - p^iy)m = (43) 

dt Tm(V) 

'J^=a,iV)il-h)-P,iV)k='^:n_J^ (44) 

dt T^h{V) 

Since we are interested in traveling solutions, it is natural to seek solutions of type V{x,t) = U{x — ct) = U{£^), 
where c is the propagation speed. To avoid boundary conditions problems, one may assume that the neuron is infinite. 
Moreover the neuron is at rest at infinity, namely we are looking for solutions such that : 

lim U{0 = (45) 

Jiboo 

The variable change = x — ct allows us to convert the partial differential equation above in an ordinary differential 
equation where ^ plays the role of a formal time: 

= -""^^^ + [QNam^'hiU - VNa) + Qru^U - Vk) + 9l{U - Vl)] (46) 
dn ,^,s,^ N ^ /, ,N n°°{U)—n 

- = M^^)(l-n)-Wn=-l^ (47) 

= a^iU)il-in)-PMn.= '^^^^^ (48) 

§ = a^m^ -h)- = ^^-^J^ (49) 

0? Th[U) 
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where we assumed for simplicity that 27rr(x) = l,Vx. 

Instead of solving these equations we shall study the corresponding equation for the FitzHugh-Nagumo model. 
They are indeed simpler and they allow us to figure out why traveling wave with a determined speed c are selected. 
The equivalent of the equations (46,47,48,49) for the FitzHugh-Nagumo model (28) are: 

e'^v + ecv + f{v,w) = (50) 
cw + g{v,w) = (51) 

where /(?;, w) has a cubic shape (e.g. /(?;, w) = v — v"^ — w) and g{v, w) is linear (e.g. g{v, w) = {v — a — bw)). More 
specifically we shall assume that we are in the situation of the Fig. 10 where only one fixed point exists for the model 
(28). In eq. (50,51) we forgot Cm and R which play no relevant role in the mechanism described below. Since ^ 

, ,2 

plays the role of a formal time we used the notation ^ = u, = il Note that the variable v, representing the local 
membrane potential, is spatially coupled by the diffusion term, while w, representing a slow ionic current or gating 
variable, is not. 

We describe the spike propagation by using the singular perturbation theory. If we set e = in the equations 
(50,51) we obtain the system of equations (called "outer equations" (101)): 



f{v,w) = (52) 
cw + g{v, w) = (53) 

As in section II. B. 2 the solution of (52) is the v nuUcline and v depends parametrically on w. The trajectory moves 
slowly on the stable branch (resp.A^~ ) and this motion corresponds to the excited phase (resp. recovery phase) 
of the pulse (see Fig. 24). 

The pulse appears then as a trajectory connecting the two branches. To characterize the dynamics between the 
two branches, it is convenient to rescale the variable ^ as | and to write (50,51) in the following form: 

dV 

V = —cv — — (54) 
ov 

w = — (v — a — bw) (55) 
c 

where we have introduced the "potential" : 



V{v, w) = — — wv (56) 

Indeed, introducing V allows us to interpret the equation (54) as the formal equivalent of the motion of a particle 
moving in a potential well with a shape V, with a "friction coefficient" c and where ^ plays the role of time. This 
picture is especially useful to understand intuitively the mechanism at work. The potential V depends parametrically 
on w and has the typical shape depicted in Fig. 21. 

When c = there is no effective dissipation and the phase portrait of the dynamical system (50) is sketched in 
Fig. 22a. In particular, there is an homoclinic trajectory connecting to itself. When c is large enough, the 
phase portrait has the shape depicted in Fig. 22b. Consequently, by continuity, there is an intermediate value of 
c, co(w) depending on w, where there is an heteroclinic orbit connecting the point and . This heteroclinic 
orbit corresponds to a moving transition layer, travelling with a speed co{w). More precisely, the heteroclinic orbit 
corresponds to an "ascending" front connecting neurons where v belongs to the — branch and with a coordinate 
^ — oo to neurons where v belongs to the + branch and with a coordinate ^ +oo (see the Fig. 23). Note that 
for each w there is a unique such cq: this is the dissipation rate required to reach asymptotically the lower bump of 
V {V~ in the case m; > 0) with an orbit starting from an arbitrary small neighbourhood of the higher bump (y+ in 
the case u> > 0) with a zero initial speed. 

Obviously the same argument can be done when w is negative. One obtains then a descending front connecting 
connecting neurons where v belongs to the + branch and with a coordinate ^ — oo to neurons where v belongs to 
the — branch and with a coordinate ^ +oo. 
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FIG. 21 Potential V of eq. (50) for : Fig. 21a : m < 0; Fig. 21h : w = 0;Fig. 21c : w > 0; 




FIG. 22 Phase portrait of eq. (54) for : Fig. 22a : c = 0; Fig. 22b : c > ;Fig. 22c : c = cq. The situation corresponds to 
w > 0. 



The complete picture is the following^. In most space the outer equations (52,53) are satisfied. When a transition 
between the two branches occurs, there is a sharp transition in v, travcUing at a speed c(i(;, e) connecting the two 
branches (and w is essentially a constant during the transition). This corresponds to a travelling pulse consisting in 
an excitation front followed by a recovery back (see Fig. 24) . Note however that the medium needs to be sufficiently 

excitable to maintain a propagation. This corresponds to the mathematical condition: ^y_'^^^}^ f{v,w^)dv > 
ensuring that there is a positive speed of propagation.. 

This picture has therefore allowed us to understand the mechanism of spike propagation in neurons, by using simple 
dynamical systems arguments. It is important to note the role of the refractory period. If the action potential reaches 
a given point, the neighboring points that have not been yet reached by the spike are depolarized to the threshold, 
while the neighboring points that have just been reached by the spike are in the refractory period and cannot emit a 
new spike. This imposes a propagation direction. 

Finally, note that the existence of travelling spike in the Hodgkin-Huxley model can also be shown rigorously 
(41; 83) For the typical values for squid axon one finds a speed value c — 21mm/ms very close to the experimental 
value found by Hodgkin and Huxley (21. 2mm/ms). 



Strictly speaking, one has still to show that this picture, obtained for e = 0, persists when e > 0. One can indeed show that the 
heteroclinic orbit persists by using perturbation theory and Fredholm arguments. 
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FIG. 23 Front corresponding to the heteroclinic connection represented in Fig. 22. 




FIG. 24 Schematic sketch of spike propagation in the spatially extended Fitzhugh-Nagumo model. 



III. NEURAL COUPLING. 

Up to now we have only considered the behavior of individuals neurons described more or less accurately by a 
set of differential equations. But neurons are not isolated entities and it is absolutely clear that the brain functions 
are the result of collective effects. If formal Neural Networks are (more or less rough) models for the brain, the 
emergent collective dynamics resulting from the coupling of individual (formal) neurons should exhibit properties 
such as information storage, recognition tasks, learning, that a lone neuron should not able to perform. If we stay at 
the level of mathematical models, then dynamical systems theory should be able to provide us some hints about the 
collective evolution when parameters are varied, external inputs are presented, learning is performed, etc .... This 
aspect are further addressed in the next sections. 

A. Synapses and synaptic plasticity 

The main function of neurons is to propagate informations via electric signals. This is reflected in their structure. 
They have two types of specific extensions: dendrites and axons. The dendrites form a tree like structure. They 
collect signals coming from other neurons and transmit them to the neural cell nucleus. The axon transmit spikes 
towards other neurons via connections called "synapses" (from Greek " syn " (together) et " haptein " (join)). There 
exists two type of synapses: electrical and chemical. In the first case (electric synapses) neurons are touching and the 
neural flux can directly go from one neuron to the other. In the second case (chemical synapses), the neurons are 
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not touching and the neural flux is transmitted vi neurotransmitters (Acetylchohn, Doparnin, Gamma- Aminobutyric 
Acid, Glutamat etc.). The action potential opens ion channels producing an influx of Ca^+, leading to the release 
of a neurotransmitter into the synaptic cleft. The transmitter difi^uses then to the other side of the cleft and binds to 
receptors, causing ion-conducting channels to open. This results in a excitatory or inhibitory post synaptic current, 
depending on the nature of the ion flow. Most synapses are chemical. 

When two neurons are connected via synapses the emission of spikes from the pre-synaptic neurons may evoke 
spikes in the post-synaptic neuron. These spikes have a variable height depending on the synaptic efficiency. Synap- 
tic efficiency evolves with time via different mechanisms. Long Term Potentiation (LTP) is a synaptic reinforcement 
mechanism involved in memory. It corresponds to an increase in the post-synaptic response after an intensive presy- 
naptic excitation, applied on a short time scale Is), but with a high frequency (> 100 Hz), inducing a strong 
depolarisation in the post synaptic neuron. Long Term Depression (LTD) is complementary to LTP. This mechanism 
arises when the pre-synaptic neuron has a low frequency activity (1-5 Hz) but the post-synaptic neuron essentially 
does not fire. This lack of synchrony between the two neurons has the effect of reducing the synaptic efficiency. It is 
believed that LTD is used in structures such as hippocampus, to bring back to a normal level of efficiency synapses 
whose efficiency has increased via LTP, rendering them available for new informations storage. A last mechanism, 
called Spike Timing Dependent Plasticity (STDP) has recently attracted much efforts. One can experimentally show 
that LTP and LTD can be elicited by carefully adjusting the timing of the pre- and post- synaptic activity. If the 
post-synaptic spike fires just before the pre-synaptic cell then the association between the two neurons weakens. On 
the opposite this association is reinforced if the post-synaptic spike fires just after the pre-synaptic cell. Important 
references for STDP studies were published in (22; 66). However, there seem to be a wide variety of different rules 
which may have different functionalities for dynamical neural networks. (3) 

B. Modeling neural networks. 

Synapses are complex objects, as neurons are. However, the more accurate one desires to model the evolution of 
a neural assembly, the less it is possible to handle analytically the dynamics. Consequently, one has to simplify the 
neurons and/or synapses description in order to obtain tractable models. Therefore, in many models synapses are 
roughly represented by a "wire" connecting the pre- and post-synaptic neuron and weighted by a number Jij modeling 
the efficiency of the synaptic connection from neuron j to neuron i. This number can be positive (excitatory synapse) 
or negative (inhibitory synapse). It can be random or constant, and may evolve in time (via learning for example, see 
sections III.C and VI. E). Although the synapses are asymmetric in general (the inffuence of j on i is not the same 
as the influence of i on j), some models consider symmetric synapses (sections III.C, V.B). Indeed, the symmetry in 
the interactions lead, for some models, to convergence properties, useful for performing tasks (see section V.B). 

Obviously, representing the synaptic conncc;tions between two neurons by an edge between two nodes is certainly a 
very rough way of sketching a neural network structure. Nevertheless, it is widely used in this community. We would 
however like to point out the following remark. Since synapses are used to transmit neural fluxes (spikes) from a 
neuron to another one, the existence of synapses between a neuron (A) and another one (B) is implicitly attached to 
a notion of "influence" or causal and directed action^". However, a neural network is a highly dynamical object and 
its behavior is the result of complex interplays between the neurons dynamics and the synaptic network structure. 
Moreover, the neuron B receives usually synapses from many other neurons, each them being "influenced" by many 
other neurons, possibly acting on A, etc... Thus the actual "influence" or action of A on B has to be considered 
dynamically and in a global sense, by considering A and B not as isolated objects, but, instead, as entities embedded 
in a system with a complex interwoven dynamical evolution. Thus the mere analysis of the synaptic graph topology 
is in general not sufficient to handle the neural dynamics. A prominent example of this is given in the section VI. F. 

On mathcnnatical gromids this aspect can be addressed as follows. Assume that the coupled neurons evolution is 
described by a dynamical system: 

du ' 

-^=Fi{uu...,UN;T) (57) 

where it, is a variable describing the "state" of neuron i (e.g. its membrane potential). N is the total number of 
neurons. P is a set of parameters accounting for neurons characteristics, external stimuli, and also including synaptic 



Note that the notion of influence roughly sketched here is very close to the definition of synaptic weights discussed by Hebb in (84) . 
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couplings (more specific examples will be given throughout this paper). In the sequel we shall use the notation u for 
the vector {ui} 

Assume now that we weakly modify the state of neuron j, for example by adding an external stimulus, such that the 
new neuron state at time t is Uj (t) + Sj (t) . The change induced on neuron i at time t + dt can be formally computed 
by writing a Taylor expansion of Fi in powers of Sj{t). At the lowest order the change will be proportional to the 
Jacobian matrix element (u) . This element measures in some sense the linear "influence" of the neuron j on the 
neuron i, when the system is in the state u. More precisely, it characterizes, to the first order in a Taylor expansion, 
the modification induced on Ui when Uj has a small variation. 

Although (57) is generally a non linear system, this Jacobian matrix can provide useful insight in the dynamical 
properties as discussed in the sections V.G and VI. F. It is in particular possible to construct a graph from the 
Jacobian matrix such that there is an oriented edge j — > i iff f^(u) ^ 0. The edge is positive if §^(u) > and 

negative if ^^(u) < 0. (Obviously, this graph depends in general on the state u). This graph has circuits or feedback 

loops If e is an edge denote by o(e) the origin of the edge and t{e) its end. Then a circuit is a sequence of edges 
ei, . . . , efc such that o(ei+i) = t(e,), Vz = 1 . . . fc — 1, and t{ek) = o(ei). A circuit is positive (negative) if the product 
of its edges is positive (negative) . A positive feedback loop basically induces (to the linear order) a positive feedback 
inducing an increase in the activity of the neurons in this loop. Obviously, there is no exponential increase since 
rapidly non linear terms will saturate this effect. 

The graph induced by the Jacobian matrix is usually distinct from the synaptic graph. In particular, it depends on 
the state u of the set of neurons. However, in models such as the recurrent neural networks discussed in the section 
V.B and VI ^^(u) is proportional to Jij with a positive (u dependent) coefficient. Thus this graph preserves the 
excitation/inhibition nature of the synapse. Nevertheless, even in this case, the mere fact that the graph of linear 
influence depends on the state of the system may have dramatic effects e.g. on signal propagation. As discussed in 
section VI. F, the notion of linear influence (and more generally linear response) allows to handle to some extent the 
interplay between the network topology and neurons dynamics and rather unexpected effects will be exhibited. 

C. Synaptic plasticity and learning. 

Synaptic plasticity occurs at many levels of organization and time scales in the brain. It alters excitability of the 
brain and regulates behavioural states (e.g. transition between sleep and wakeful activity). It is also involved in short 
and long term memory and learning. In this section and in this paper we shall only focus on this last issue. 

The synaptic weights are evolving in time during learning. In formal neural network learning is thus represented 
by evolution schemes for the synapses, called learning rules. Although learning rules can be proposed using precise 
description of LTD, LTP and STDP, most of them rely on some fundamental recipes inspired from D. Hebb's work. 
One speaks then of Hebbian learning. We shall focus on Hebbian learning in this paper. 

D. Hebb has proposed in (84) a theory of behavior based on the physiology of the nervous system. The most impor- 
tant concept to emerge from Hebb's work was his formal statement (known as Hebb's rule) of how learning could occur. 

When an axon of cell A is near enough to excite a cell B and repeatedly or persistently takes part in firing it, some 
growth process or metabolic change takes place in one or both cells such that A 's efficiency, as one of the cells firing 
B, is increased. 

Most of the learning rules in neural networks are based on Hebb's observations plus a few well established facts. 
They rely upon a few recipes that can summarized as (93): 

• Learning results from modifying synaptic connections between neurons. 

• Learning is local i.e. the synaptic modification depends only upon the pre- and post- synaptic neurons activity 
and does not depend upon the activity of the other neurons. 

• The modification of synapses is slow compared with characteristic times of neuron dynamics. 

• If either pre- or post- synaptic neurons or both are silent then no synaptic change takes place except for 
(exponential) decay which corresponds to forgetting. 

The first item implies that learning results in a modification of the J^ 's. The second one basically says that the 
synaptic modification of Jij writes J^j = eh{J^, mj,mi) where J^^ is the value of the synapses j ^ i after the learning 
rule has been applied. The parameter e has been added for convenience and will be discussed below. The numbers 
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rrii (rrij) denotes the "state" or "activity" of the neuron i (j). We do not precise yet what is this "state" since it can 
vary according to the modeL Several examples will be discussed below. The third item implies then that e is small 
parameter, whose inverse corresponds to the characteristic time for a significant change of Jy . The fourth item may 
lead to different forms according to the model (see below). But if one assumes that the changes in the Jjj's are slow 
(item 3) and if /i is a smooth function then one may simply consider a Taylor expansion of a generic regular function 
h. This gives, up to the second order in m^, mj. 

J'ij = £ (aooo + aimJij + aowrrij + aooirrii + aaumimj + h.o.t.) 

where h.o.t. means "higher order terms" such as JijinirUj, etc.... In this chapter we shall focus on this form, 
forgetting the other terms. Note that the terms oioo, aoiO) aooi) aoii have all a "biological" interpretation. We 
shall not consider the term aooo • Writing A = eaioo the corresponding term models passive "forgetting" : if a 
synapse is not solicited its intensity decreases with a decay rate j (we shall assume that 1 > A > 0). On biological 
grounds, the situation is a little bit more complicated. The decay of the synapse and more generally its evolution 
depend on the activity of the pre synaptic (j) and post synaptic (i) neuron, as we saw. These activities determines 
the production of Ca?^ ions, which acts in turn on the width of ionic channels involved in the synapse activity. 
The production of Ca^"*" increases whenever i and j are "active" increasing the synaptic efficiency. On the other 
hand, when Xi or xj are active then the concentration [Ca^+] stays constant, and enzymatic phenomena result 
in an effective decay of the synapse (Long Term Depression) . This gives an interpretation of the 3 terms aoio , aooi , aoii ■ 

Thus, setting eaioo = A, eaon = a, eaoio = eaooi = ~7 we obtain a synaptic evolution having the form: 

Jlj = XJij + aTij — Prrii — jrrij (58) 

rfj is a function of the activity of the pre- and post- synaptic neurons. In most case Tfj ~ mfmj but the form (III. A) 
affords natural generalization that we shall briefly discuss. Note that all the coefficients a, f3, 7, A are proportional to 
e, which flx somehow the characteristic time scale of the synaptic dynamics. 

Some examples of learning rules will be presented in this chapter but we shall focus on situations where /? = 7 = 0. 
A more detailed discussion can be found in chapter III. 

IV. WEAKLY CONNECTED NEURONS. 

What happens when neurons, having their own dynamics, are coupled via synapses ? Though this question is 
too general to have a precise answer, it is possible to address it when considering a weak coupling limit with some 
additional assumptions discussed below. In a nutshell, the basic idea is to consider the situation where a collection of 
neurons is coupled as a perturbation of the uncoupled case, where each neuron evolve independently from the other. 
The perturbation resulting from the coupling can however be cither irrelevant, when the coupled and the uncoupled 
systems are essentially equivalent from the dynamical point of view (section IV. B), or it can have a drastic effect. As 
argued below, this is basically the case when some neurons are close to a bifurcation point. In this case a rather detailed 
analysis can be made by using standard tools from bifurcations theory theory such as center manifold reduction and 
normal forms (sections IV. C and IV. D) provided one restricts the overwhelming possibilities of bifurcations, that may 
potentially occur in a collection of coupled neurons, to some canonical "scenarios" (sections IV. E, IV. F). This is, of 
course, an important restriction, but the results obtained are quite illuminating from many aspects, especially with 
respect to the ability of such Neural Networks to perform task, such as pattern recognition, manifested by changes in 
the dynamics when a pattern is presented to the network (section IV.G). We present here a short review of results 
mainly due to Hoppensteadt and Izhikevich (see (93)). 

A. General setting. 

From now on, we consider therefore an assembly TZ oi N coupled neurons. Summarizing the previous section, the 
dynamics of individual neuron is governed by an equation of the form: 



^=F,(X,;A), z = l...N 



(59) 
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where Xi is a vector in R™ describing the state of the neuron. Tliis form is quite general and includes in particular the 
Hodgkin-Huxley equations (8-11) and the general form of excitable membrane equations (17). A is a set of parameters 
on which the neurons dynamics depends. An example is the applied current I. We may assume that A belongs to a 
space £\ . Without loss of generality, and for technical reasons we shall assume from now on that Xi Q Ai where M. 
is a compact m dimensional manifold. 

The basic requirement of the theory of weakly connected neural networks is that the contribution of activity of one 
neuron to the activity of another one is very small. More precisely, following (93) we call Weakly Connected Neural 
Network (WCNN) a dynamical system of the form: 

^ = Fi{Xi;X) + eGi{Xu...,Xn;X,p,e); Xi e M, i = l...N (60) 
or, in a more compact form: 



— =F(X;A) + eG(X;A,p,e) (61) 

where we note in boldface the n N x m dimensional vectors X = {Xi)f^^ , F = (Fi)^^ , G — (Gi)^-^. Note 
therefore that F as a diagonal structure. In equation (60) Gi is a smooth, bounded function of {Xi, . . . , X^; A, p, e) 
that models the synaptic connections between the other neurons Xi, . . . ,Xj\j and the neuron i. It depends on the 
set of parameters X € £x C describing the state of individual neurons, on the coupling parameter e, and on an 
additional set of parameters p G £p C corresponding to external constraints (for example the external environment 
influence, a static input, etc . . . ). Finally, it is assumed that e is small, namely e <C 1. This purely "mathematical 
assumption" is required to perform the analysis presented below. 

Of course, one may wonder whether this restriction still provides models for biologically realistic situations. Note 
however that the mathematical condition e ^ 1 is abstract and refers to the particular set of differential equations 
(60) which attempts to model some aspects of neuronal dynamics. Henceforth, questions such as : "How small is e 
in the real brain" are essentially meaningless. There are nevertheless different and non equivalent ways to estimate 
the strength of connexions between neurons, One of them is based on the analysis of cross correllograms from pairs 
of neurons. Performing such an analysis Abeles (1) concluded that interactions between adjacent neurons in the 
cortex are weak and the interactions between distant neurons is even weaker. Another way to characterize weakness 
of synaptic connections is to measure the amplitude of post synaptic potentials (PSP) in the soma while the neuron 
membrane is far below the threshold value. Indeed, in this state the size of PSP reflects the weakness of connections. 
A detailed discussion can be found in (93). 



B. Structurally stable case. 

In spite of the restriction to weak couplings, the dynamics of (60) can be very rich since, as we have seen in the 
previous sections, the behavior of individuals neurons can already be quite complex. Consequently, it is impossible 
to analyze (60) without further restrictions or specifications. A starting point is to consider first the case when each 
neuron has an equilibrium point and is in a rest state when e = 0. In order to simplify the computations, and without 
loss of generality, one may assume that this point is the origin when A = 0, namely F(0,0) = 0. Denote by: 



L =^£)xF(0,0) = 



/ii 

L2 



V 



\ 





(62) 



J 



the Jacobian matrix of F at the point X = 0, A = 0, where Li =^ (§Y^j ■ A fixed point may be hyperbolic 

or not. The following result derived by E. Izhikevich (96) is an example of application of the Hartmann-Grobman 
theorem (see appendix) in the context of Neural Networks dynamics. 



Theorem 1 // the dynamics of each neuron is near an hyperbolic equilibrium when e 
the weakly coupled one (61) and the linear system : 



then the uncoupled network, 
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^ = LX (63) 
dt ^ ' 



are topologically conjugated. 



This means that the entire coupled neural network is essentially a linear system and is not more complex. If we 
admit that non linearity plays a fundamental role in neural dynamics, the corollary leads to the following conclusion. 
A WCNN with an hyperbolic fixed point is not interesting as a "brain" model The corresponding dynamical system 
being structurally stable, wc arc led to the bewildering conclusion that a neural network model where each neuron 
has an equilibrium point and is in a rest state when e = needs (at least) to be structurally unstable in order to 
exhibit (real) non-linear effects and to be relevant for brain dynamics. 

At this point we may pose. Indeed, it might be difficult to find out a biologic;ally realistic situation where the 
neural network is in a state corresponding to a fixed point. There is at least an evident one : death (this is indeed a 
structurally stable situation). But this is not a very interesting example. This led us to several remarks. First, as 
shown below, even if the uncoupled system consists of neurons in a rest state, the coupled dynamics can be quite a 
bit more complex, even if the coupling are weak (see section VI. C). As soon as neurons are coupled, many different 
situations and dynamical regimes may occur. We shall see several examples in this chapter, from periodic to chaotic 
regimes, with one or several attractors etc... Consequently, coupling neurons with a rest state corresponding to a 
fixed point docs not necessarily mean that the coupled dynamics will be at rest. 

A more important issue concerns hyperbolicity. Though the theorem 1 deals with hyperbolic fixed points, the notion 
of hyperbolicity extends to quite more general attractors than fixed point, such as strange attractors (see appendix). 
If we are interested in the ability of a Neural Network model to perform tasks such as recognition of an external 
pattern, and if we agree that this recognition corresponds to some dramatic change in the dynamics, then one can, 
in principle, extend the wisdom coming from theorem 1: to be dynamically reactive to solicitations from the outside 
the neural network needs to be close to a point in the parameter space where the dynamics is structurally unstable 
with respect to perturbations corresponding e.g. to a specific (learned) pattern. To be efficient and adaptable the 
system needs to be close to a critical point in order to display punctuated response to external world changes. We 
shall actually propose an example in section VI. E exhibiting a behavior that can be related to this statement. It 
might well be that this conclusion extends more generally to biological systems and to living systems (see (19), (20)). 



C. Central manifold reduction. 

As wc shall sec right now, the situation is already quite a bit richer when some of the neurons have a rest state 
(corresponding to a fixed point) which is not hyperbolic. Assume therefore that there is a subset B GTZ oi neurons 
such that the Jacobian matrix Li of the uncoupled neuron i, i G B has eigenvalues with a zero real part. We shall 
call these neurons critical since they are close to a bifurcation point. This means that the slightest change in the 
parameters, induced either by the inputs of other neurons, or by an external stimulus, etc . . . , may provoke a non 
linear dynamical response of the neuron such as spike, train of spikes, etc... Moreover we order the neurons such that 
the k first neurons are critical. 

Even when focusing on this situation, the analytical study of the changes occurring when the coupling e is switched 
on, and when the parameters p. A are modified is not tractable in general. Indeed, in the most general situation, each 
matrix Lj have a number rij < m of neutral eigenvalues, some neurons may be at the threshold for a Hopf bifurcation, 

some others close to a saddle-node bifurcation, etc Moreover, if m, the number of control parameters, is sufSciently 

large (depending on the accuracy of the underlying model) one may have bifurcations of codimension larger than one 
or two. It is therefore natural to start from the simplest situations, namely the case where all neurons undergo 
the same (codimension one) bifurcation. Therefore, we focus now on the case where either the matrix Li of each 
critical neuron has a simple zero eigenvalue, or a pair of complex conjugate imaginary eigenvalues. In these cases the 
techniques of central manifold projection and normal form reduction (see appendix) allows us to reduce the dynamical 
system (60) to a canonical form, close to the bifurcation point, and provided e is sufficiently small. 

The result presented below is rather abstract (though it is a direct application of the center manifold theory (30) 
adapted to the present context). But it has interesting consequences discussed in the next subsection. Basically, this 
results shows that the dynamics of the coupled system is locally governed by the critical neurons. 



More modestly, one may consider, instead of the brain, small functional units such as cortical columns or simple nervous systems 
(worms). Fortunately, the same conclusion holds. 
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Theorem 2 (Izhikevich (96)) Suppose that each of the first k Jacobian matrices Li, i — l...k is non hyperbolic. 
Then the dynamics of (60) is locally governed by a dynamical system of the form : 

^ = fi{xi\X) + egi{x]\,p,e), i=l...k (64) 

where}"^ Xi G E], and Ji = Dxji{z, 0) = L^^e-; ■ Moreover, if = {0}, there is a function Z : i?'^ x 5;^ x 5^ x R — > Ai"' 
such that any local solution X(t) of (60) close io X = tends exponentially to Z(x(t), A, p, e) where x(t) is a solution 
of (64). 



D. General normal form 

Once the center manifold reduction has been done and once one has identified the type of instabihty occurring for the 
uncoupled neurons one can further reduce the dynamical equation (64) to a canonical form or norm,al form, allowing 
somehow to classify the models into equivalence classes. For this, one needs however to provide some additional 
restrictions on the type of couplings (function Gj), and on the parameter dependence (transversality conditions). 

In the more general case, e, p, A are independent parameters and the situation is quite complex. A way to simplify 
it is to assume that A, p have the following form: 



A ~ A(e) = + eAi + e^M + 0{e^) ; Ai, A2 e (65) 
P ~ P(e) = Po + epi + o(e^) ■,po,pie£p (66) 

such that A(0) = and p(0) = po- This form is convenient to handle but it is not a loss of generality since Ai, A2 (resp. 
Po,pi) are still independent control parameters and can assume any value in £x (resp. £p). As discussed above, A 
corresponds to control parameters allowing for example to tune the neuron characteristics, while p mimics "external 
constraints". In the form (66) po may correspond to some "background" influence po while pi models an external 
input^^. 

For e = the uncoupled neurons have a control parameter A = 0, the critical neurons are located at the bifurcation 

point, and the influence of the external environment is modelled by the parameter po- When the coupling is switched 
on, the physiological parameter A of individual neurons is modifled unless Ai = 0. Similarly, the external world 
influence is manifested by an additional term epi superimposed upon the "background" influence po. A particularly 
interesting issue is to study the reactivity of the set of coupled neurons with respect to an external input, modeled 
by the parameter pi . If one thinks for example of recognition task one expects a non trivial sensitivity to the input 
pattern pi, manifested by qualitative dynamical changes. E. Izhikevich analyzed this situation in great details and 
made a classification of the normal form for the coupled system according to the type of bifurcations the individual 
neurons are close to. 

Let us give the main ideas and results. Assume that we are in the situation of the theorem 2 and that the center 
manifold reduction (64) has been performed. Using the form (65) (66) for A, p one rewrites the dynamical system (64) 
in the form : 

^=f^ eAi + e^Aa + ©(e^)] + egi [x; eAi + e^As + 0{€^), po + epi + o{€^), e] = 
= fi{xi; 0) + e [Dxfiixi, 0).Ai + 5i(x; 0, po, 0)] + 



d 

Dlfi{xi,0).{Xi, Xi) + Dxfi{xi,0).X2 + Dxgi{yi;0, po,0).Xi +-Dpfifi(x;0,po,0).pi + —gi 



The tangent space of each Jacobian matrix Li can be decomposed as: 

iR™ = Ef e e Ef = e Ef 

where E"^ , E'r arc respectively the stable, unstable and central space. 

A similar description is made in the model of section VI. The microscopic parameters of the model (102) corresponds to A while 
the external input ^ of section VI. E corresponds to pi. Note however that the analysis performed in section VI does not require the 
assumption of weak coupling and closeness to a rest state. 



33 



Close to the fixed point x = tliis reduces to: 



= ^ (0; 0)xi + hixf +eai + e^di + e2_^ SijXj + o(e^) (67) 



where : 



ai = Dxfi{0;0).Xi+gi{0;0,po,0) (68) 
1 ^2 f,- 

hi = 2 9^(0' 0) (69) 

di = i?^/i(0;0).(Ai,Ai)+DA/i(0;0).A2 + £'A5i(0;0,po,0).Ai+£'p5i(0;0,po,0).pi + ^5i (70) 

Sij = |^(0;0,po,0) (71) 

oxj 

It is particularly remarkable that the coefficient Sij acts as a "formal synapse" coupling the neuron j to the neuron 
i. But, contrarily to the usual synapses, which establish a "wired" link between two neurons, the coefficient Sij 
is essentially generated by the (nonlinear) dynamics. Note that it is given by the Jacobian matrix of g. This is 
thus a first illustration of the concept introduced in the section III.B. It corresponds to an effective link that is not 
necessarily supported by a wired link. This fundamental aspect is discussed below and in more details in the section 
VI. F (see in particular the equation (128)). 

Up to now we have written general equations without consideration about the (codimcnsion one) bifurcations of 
the critical neurons. To each type of bifurcations is associated a set of transversality conditions (see appendix) that 
allows to reduce further the equations (67). 

The result (67) is rather abstract but it has interesting development in particular in the following case. Assume 
that we slightly perturb the dynamical system with an external input corresponding to the parameter pi. What is 
the effect of this perturbation on the dynamics ? In particular, assuming that some neurons are close to a bifurcation 
point, does a perturbation of the form (65) have an effect on the global dynamics ? Such a change can be considered 
as an effective reaction of the system to the external input, that can be used, for example, to perform recognition 
tasks. Consequently, it can be useful to have analytical results on the normal form of the coupled dynamics near the 
bifurcation point. We shall not discuss all the cases investigated by Hoppensteadt and Izhikevich. We shall focus 
instead on two examples that we find particularly rich and enlightening. 

E. Saddle-node and pitchfork bifurcations. 

Let us first analyze the case when the k critical neurons are close to a saddle-node bifurcation. Then, for each 
critical neuron i = l...k the uncoupled vector field fi{xi,X) satisfies the transversality conditions (see appendix) 
= 0; |^/<(0;0) 7^ and _Da,A(0;0) 7^ (which means that Da,A(0;0) has no zero eigenvalue). Therefore, it 
follows from (67) that the normal form of a WCNN where k neurons undergo a saddle-node bifurcation is: 

^ = hixf + eai + e^di + e ^ SijXj + o(e^) (72) 

i=i 

The variable change Xi e^Xi and the time change t e~^t transforms these equations into: 

^ =ai + hixl + ^^^SijXj + edi + o(e^) (73) 
i=i 

We now want to analyze the effect of an external input pi on (61). From eq. (68) ai depends only on pQ, A while 
hi is independent on these parameters. Moreover, if a, ^ 0, the dynamical system (73) admits, to the zeroth order 

in e, the fixed points ±Jj^, which are hath hyperbolic. Since hyperbolicity is structurally stable, when ai ^ the 
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dynamical system (73) is insensitive to the input pattern pi (and, at least in this sense, cannot perform recognition 
task). It can react only to po and this reaction is trivial (by the implicit function theorem). 

Consider now the case when a, = 0. Then the variable change Xi = eh^^yi and the time change t = et in (72) lead 
to: 



^=ri + x^+f;^CijXj+0{e) (74) 



where: 



n = hidi (75) 
Cij = hiSijhJ^ (76) 

The dynamical system (74) is the normal form of (73) under the condition ai = 0. It depends now on the external 
input pi via dj. The fixed points of ((74) are now determined by : 



U + Xi 



(77) 



Consequently, they depend on the structure of the matrix C = {cjj}"j_i (corresponding to the effective links induced 
by the nonlinear dynamics). The condition = writes : 



DxMO; 0).Ai + 3,(0; 0, po, 0) = 0; Vi = 1 . . . n. (78) 

This specific condition is called by E. Izhikevich, the adaptation condition. He showed in particular this interesting 
result: in order for the reduced dynamical system (61) to exhibit a non trivial sensitivity to the input pattern pi, one 
needs that the neurons adapt to the pattern, via the internal parameter A. In other words, A is not independent on p 
but it must satisfy (78) which means that the internal parameter A counterbalances (up to order e) the steady state 
input from the entire network onto each neurons. This notion has deep implications, because it suggests that a suitable 
training of a Neural Network such as (61) requires an evolution of the control parameters X, under the influence of 
the external input pi, with the constraint that the condition (78) is achieved. In some sense, learning is efficient if 
the training leads the system into a very specific part of the parameters space, corresponding to the condition (78), 
where the system is close to a bifurcation point and where this bifurcation is only induced by the presentation of the 
learned pattern (or a weakly perturbed version of it). This has interesting echoes with the discussion of the effect of 
hebbian learning in the model described in the section VI. E. 

The adaptability condition is also important if one considers other types of bifurcations such has the pitchfork (in 
this case the transversality conditions impose that (78) is automatically satisfied) and the cusp bifurcation. 

In the case of pitchfork bifurcations, for example then the normal form is : 



^ = bxi±xl + J2(^ij^j (79) 

i=i 

We shall return back to this form in the section IV. G where we shall consider the effect of having a synaptic matrix 
constructed via Hebbian learning. 



F. Hopf bifurcations. 

Assume now that the uncoupled system has k neurons close to a Hopf bifurcation. This means that the Jacobian 
matrix of each corresponding critical has a pair of purely imaginary eigenvalues ±0j at the critical value of the 
parameters. Call Vi,Vi (resp. Wi,il]i) the corresponding right (resp. left) eigenvectors. Moreover, the vector field 
satisfies the transversality conditions THl in the appendix. Using similar techniques as in the previous example one 
can prove the following (93): 
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Theorem 3 // the dynamical system (64) is near a multiple Hopf bifurcation then there is a variable change and a 
time rescaling t = et reducing it to the normal form: 

dz- ^ 

= hiZi + diZi\Zi\^ + dj Zj + 0{y/^) (80) 

where bi,di are complex coefficients and where the coefficients Cij are given by : 

f~, _ { w^.D^.gi.Vj ifn, = , . 

tfn.^Qj ' 

This results is quite interesting. It shows that close to the bifurcation point the oscillating neurons can be divided 
into pools according to their natural frequency. There is an effective coupling Cy between the neurons in the same 
pool, while the coupling between neuron from different pools is negligible (namely weaker than 0(\/e)) even if there 
is a synaptic connexion between them in the global system. This strongly suggests that resonances appear in such 
system. Indeed, one can exhibit Arnold like tongues in this situation (see (93) page 173). Also, a periodic input with 
a frequency ut one can establish interactions between oscillators transmitting on different frequencies or it can disrupt 
interactions between oscillators transmitting on the same frequency. 

This result as strong implications going far beyond the field of Neural Networks. Indeed, it suggests that non linearity 
may induce effective paths among units that are not directly connected to the graph of interactions. Moreover this 
type of behavior is not specific to weakly coupled neural networks with a rest state close to a Hopf bifurcation. In the 
section VI. F we shall exhibit a similar behavior for a chaotic Neural Network and we shaU show that a linear response 
theory based on recent results by D. Ruelle (131) might be used to locate these resonances. 

G. An example of "Hebbian" learning. 

In this section we give a first example of Hebbian learning rule allowing the neural network to perform tasks such 
as pattern recognition. We start from the "generic" form of Hebb's rule (58) with /3 = 7 = and where the "state" 
of the neuron i (m, in eq. (58)) is given by Xi and where Tij = XiXj. 

J'ij = XJij + aXiXj + higher order terms. (82) 

Suppose now that we have several "patterns" or "images" to be "memorized" by our neural network. 

What means "memorized" ? Many definitions are possible but, in the context of this chapter, we shall consider that a 
pattern is memorized if the neural networks has acquired, via learning, the capacity to dynamically evolve towards a 
"state" "associated to the pattern" , provided that it was "suitably prepared" . We insist on the fact that this property 
must be acquired, namely, it should not exist when no learning is performed. This definition is however still rather 
ambiguous. What is a "state" ? What means "can be associated" to the pattern, "suitably prepared" ? Again, there 
are many possible interpretations. But let us start with a very simple one. Assume that the "patterns" are vectors 
corresponding to points in the phase space of our neural network. Assume then that, after learning, the dynamic 
evolution admits all patterns as stable fixed points. Then, starting from an initial condition in the attraction basin of 
the pattern k, say, the dynamics will converge to this pattern. Here "suitably prepared" means that we start from 
an initial condition in the attraction basin and the "state associated to the pattern" is the pattern itself and this is 
fixed point of the dynamics. Since an initial condition in the attraction basin of the pattern is a (possibly small) 
perturbation of it, one may interpret the convergence as a recognition of the pattern by the neural network when a 
perturbed version of it is presented. 

How to manage such a "learning ability" ? Consider the equation (82) and assume first that the coupled neural 
network admits only a stable fixed point corresponding to the pattern Then the Hebbian rule implies that 
Jij ^ possible generalisation to p pattern is then: 

fe=i 

We shall actually see, in the section V.B, the effect of this rule in a recurrent neural network. We shall also show 
another application of Hebbian learning in a situation where the recognition of a pattern does not correspond to the 
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convergence to a fixed point, but to a more complicated object in tlie phase space (section VI. E). For the moment 
and to stay in the spirit of the section IV we will address the following question. What is the effect of the hebbian 
synapses (83) in a neural network where some neurons are close to a bifurcation ? 

Let us for example consider the case of a pitchfork bifurcation and assume therefore that the dynamics is given by: 

^ = bxi-xf + f2Jij^3 (84) 

corresponding to having a STibsct of neurons close to a pitchfork bifurcation (see section IV. E). It is easy to see that 
X* = is always a fixed point and that its Jacobian matrix is given by -DFx* = bl + JJ, where / is the identity 
matrix and the matrix of synaptic couplings. Assume now that J' is given by eq. (83). Note that, according to 
this equation, J' can be written in the form J — Y%=\ where $, denotes the transpose of ^. Assume now that 
the patterns are orthogonal and take binary values = ±1 (thus {$,^,$}) = N6k,i)- The mutual orthogonality 
of the patterns imposes that the number p of patterns is lower than the dimension of the phase space {p < N). 
Then it is easy to see that J'^'^ = Pk^'^. Thus f3\,. . . ,(3p are eigenvalues of J with corresponding right eigenvectors 
. . . , The remaining N ~ p eigenvalues are all zero and the corresponding eigenvectors belong to the orthogonal 
of span . . . , (the kernel of J'). If we finally assume that /3i > f32 > ■ ■ ■ (3p then we see that DFq has p 
eigenvalues 6 + /3i > 6 + /32 • • • > b + Pp. Consequently, looses its stability when b > —Pi. This destabilization 
arises via a pitchfork bifurcation occurring in the direction Indeed, set = yi^,} then ^ = (6 + Pi)yi — y\. 
Hence, two symmetric fixed points Xj_ = ±-^6 + appear, proportional to the pattern If we further increase 
b we have a sequence of similar pitchfork bifurcations, in the direction ^'^j whenever b > —Pk. Finally, a stability 
analysis of these new fixed points show that if 6 > —Pm + ^^^^"^ , there are m pairs of stable nodes corresponding 
to the patterns . . . In this sense, the Hebbian rule (83) gives to a neural network of type (84) the capability 
to retrieve memorize patterns corresponding to stable fixed points. We shall return on a more general version of cq. 
(84) in the section V.B. Note finally that if one increases b beyond a positive value then "spurious" memories appear, 
that is, stable fixed points that do not correspond to any of the memorized patterns. These new patterns belong to 
KerJ (see (93) for details). 

V. RECURRENT MODELS. 

In the previous section we have considered the effects of weakly coupling neurons, in a situation where uncoupled 
neurons arc close to a bifurcation point. We have presented some rigorous results obtained from bifurcations and 
normal form theory. They reveal some illuminating aspects of the emergent dynamics, such has the adaptation 
principle or the existence of an effective network induced by the dynamics and not necessarily identical to the synaptic 
network. 

We depart now from this setting. We want to analyze the collective dynamics of a neural network where the 
couplings are not necessarily weak and where neurons are not necessarily close to a bifurcation point. For this, we 
consider a recurrent neural network whose dynamics is given by the equations (87) (continuous time) or (102) (discrete 
time). One can indeed go quite a bit deep in the dynamics description. Moreover, the model presents an overwhelming 
richness and it can be partially analytically studied. This is also a good benchmark for developing tools in non linear 
networks analysis (see section VI. F). 

We first show how this recurrent model can be derived by switching from a spiking description of the neuron to a 
frequency rate description (section V.A). We then discuss the dynamics of the model when the synaptic weights are 
symmetric (section V.B). Gcmcral results from dynamical systems theory allow one to prove a convergence property 
and to exhibit a Lyapunov function (see appendix for a definition). This function has some analogies with the magnetic 
energy in a system of interacting spins. Actually, the "energy" landscape presents a structure similar to the rich and 
complex energy landscape of spin glasses models. There exist a large number of minima (fixed points), whose number 
increases exponentially with the system size. The existence of these many minima is closely related to the competition 
excitation/inhibition induced by synapses, and corresponding to frustration in spin glasses. Actually, the techniques 
developed in statistical mechanics of spin glasses can be adapted to estimate the number of minima (111) and to have 
an good description of energy landscape. The convergence to minima can then be used to store informations if one 
uses the Hebbian rule (83) (92). 

In section V.C we present briefly cooperative systems, where one still have a convergent dynamical system even 
if the synapses are not symmetric. We discuss in particular shortly a fundamental result by Hirsch having recent 
extensions in the field of genetic networks (74), (142). But usually, when symmetry is broken, the dynamics is not 
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convergent and can be, for example, chaotic. The section VI is entirely devoted to a chaotic model, which has, 
furthermore, nice properties when submitted to an Hebbian like learning. 



A. From spiking neurons to firing rate neurons. 

In section IV we have discussed the effect of coupling neurons without giving a detailed description of the individual 

neuron dynamics. For further development we have now to specify it. For this purpose and in the spirit of table 1 
we switch now from neurons having an activity described in terms of spikes to neurons described in terms of firing rate. 

Assume therefore that each neuron is emitting spike trains and call Xj{t) the firing rate of the neuron j at time t. 
Note that the definition of firing rate requires an integration over a certain time window, (see the introduction), that 
wc shall assume to be short compared to the time scale for the evolution of the variables considered in this description. 
Call Ui{t) be the time average of the membrane potential on this time window. In standard models the firing rate is 
a function of uf. Xi{t) = fi{ui{t)) where /, is a sigmoidal function such as^^ fi{x) = Erf{giu) or 1±1e}^(9±h1_ g. [g 
called the 'gain' of the transfer function /j. Since the slope at the inflexion point of fi is proportional to Qi and since 
oo) = and /i(+oo) = 1 this parameter measures the level of non linearity of the function. The sigmoidal shape 
of fi can be understood by the following argument. 

A spike is emitted each times the average membrane potential exceeds the neuron threshold 6i. This threshold, as 
we saw, depends on time. In particular, after emission, it increases to infinity during a time Ta, assumed here to be 
identical for all neurons, corresponding to the absolute refractory period. Then it decreases to reach its initial value. 
Hence, if 6{t) is the threshold value at time r, the initial time being the time where a spike is emitted, one has: 



61(7 



00 if 
decreasing function if 



< r < Ta 
T> Ta 



(85) 



At each time Tj such that 9i{Ti) = Ui there is an emission of a spike, therefore the corresponding average time of 
firing (the initial time being the time where a spike is emitted) is given by Tj = 6^^{ui) and the (normalized) frequency 
rate is : 



< 1 



Consequently, Xi is an monotonously increasing function of in with values in [0, 1]. In the case of integrate and fire 
neurons driven by an external stochastic current I{t) one has has an explicit equation for /j. The membrane potential 
evolution is given by: 

Tm-^ = —iu(t) + RI(t) 

(cq. 36, section II. B. 4). Assume that I{t) is random (e.g. Poisson process) with a (stationary) probability distribution 
V . Then the probability that the neuron fires at a time t -\- dt \s given by 



V [u{t + dt) > 



7 , ^ . ^ RI{t)dt 

(1 - —dt)u{t) + — > 

'T'm. 'T'm. 



= V 



I{t)dt >^{e- u{t)) + u{t)^dt 



where is the repartition function of /. 



One also finds in the literature the case where fi{u) = taBh{giu). Hence / takes its values in [—1, 1]. 
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Consider now an assembly of such neurons. The neuron i receives the spikes coming from otlier neurons, and tlie 
total current Ii{t) is the sum of spikes arriving from each neuron j, weighted by J^. The membrane potential Uj of 
the neuron i depends on the frequency rates of the spikes trains emitted by the neurons connected to i. The current 
received by i is therefore X^^i JijXj{t) where xj can be viewed as an integration over a small time window of the 
current appearing in eq. (38). Then the analog to the evolution equation (36), section II. B. 4, is : 

N 

-^ = -u,{t)+J2j^,^At) (86) 

where we have dropped the time constant in front of — Uj, incorporating it into the time scale dt. 

Fixing some threshold reference (e.g. 6^ = 9{2Ta), setting Ui{t) = Ui{t) — 6^ and Xi = fi{ui) one finally obtains: 

-^=-u,{t) + J2JiofoMt))-^i: i = -^...N. (87) 

The model (87) displays a wide variety of dynamical behavior according to the form of J the matrix of synaptic 
couplings. 

The equation (87) is a particular form of the Cohen-Grossberg model (44). The general form is: 



dui 

-dt = "'^"^ 



N 



i=l...N. 



(88) 



where ai,6j are mild functions (e.g. is bounded, positive and locally Lipschitz continuous and 5^,6^ ^ are locally 
Lipschitz continuous). In the sequel we shall however restrict to the model (87). 



B. Symmetric synapses. 

Consider first the case with symmetric synapses Jij = Jji. One shows then that (87) is convergent whenever J is 
symmetric. 

Here is an elegant proof due to M. Benaim(24). It has the advantage to hold in more general case than for symmetric 
synapses. M. Benaim proved indeed the following theorem: 

Theorem 1 (24) Consider the differential system: 

flu ■ 

^=h{u)Gi{u)=Fi{u), i = l...N (89) 

where hi : — > R*"*" are strictly positive functions and assume that there exist a family of strictly positive 
functions tpiiR^ R*~^ such that the following holds ("detailed balance" condition)^^ : 

Then : 

1. (89) admits a strict Lyapunov function. 

2. The isolated equilibria of (89) are generically hyperbolic. 



The name "detailed balance" comes from the evident analogy with the equilibrium conditions for the stochastic evolutions, such as 
Glauber dynamics, used in statistical physics. 



39 



The proof rehes on the remark that the differential form lo = '}2!i=i '^i{ui)Gi{u)dui is exact {dui = 0) since: 



l<3 

It follows that there exists a function V such that w = —dV. Consider now the Riemmanian metric defined by: 

^^^y)^ = L.-^^^vi (91) 

Then dV = — < F,du >u and consequently ^ = — < F,F >„. Hence V is strictly decreasing along the 
trajectories of the dynamical system (89). Thus (89) is a gradient field for the previous metric and F is a strict 
Lyapunov function. Consequently, from the Lasalle invariance principle (see appendix) the dynamical system is 
convergent and the fixed point are generically isolated. The Jacobian matrix being self adjoint for this scalar product 
its eigenvalues are real and are generically non zero. 

In the case of the dynamical system (87) one can write: 

def Gi(u) 



dui 1 



N 



(92) 



(since vanishes only at infinity, one may assume that the initial conditions are chosen in some compact set of IR^). 
Then bi{ui) = and : 

0^ = -Sij {{ui + 6i)f"{ui) + f{ui)) + Jijf'j{uj)fl{ui) 

The detailed balance condition holds if Jy = Jji and it writes = (hence one may take -0^ = 1 in the 
scalar product (91)). Thus, (87) is a gradient system with a Lyapunov fimction V given by dV — — X^i^i Gi{u)dui = 
'^f-l{ui{t)f^{ui) — X^^i Jijfj{uj)fi{ui) + 6i)dui. This gives, up to an irrelevant additive constant: 

y (u) = ^ / {u + ei)f[{u)du - ^ Jiifj {uj)Mui) (93) 

This result has several important consequences. First it shows that (87) is convergent (44). But it gives substantially 

more. The equilibria are the minima of the function V. Actually, this function looks very much like the magnetic 
energy in physical systems. Assume indeed, that the slope of the sigmoid tends to infinity. Then / becomes a binary 
function. If f{u) = tanh{gu) then Xi = f{ui) takes value in {—1,1} as the binary spins of the Ising model. The 
Lyapunov function (93) writes, in this limit, — ^^^I'hjXiXj (since f'{v) everywhere but at m = where it 
becomes infinite). The Lyapunov function as therefore exactly the form of the energy resulting from the magnetic 
interaction between Ising "spins" . 

The structure of the energy landscape of magnetic system with ferromagnetics (J^ > 0) and antiferromagnetics 
( Jij < 0) magnetic interactions is astonishingly complex. Actually, when the Jy 's are randomly distributed (modeling 
the presence of impurities in a magnetic sample) one obtains a model for a spin glass (for a review see for example 
(26), (114), (141)). Spin glasses have extremely rich properties and the canonical models (such as the Edwards- Anderson 
model (61) or the Sherrington-Kirckpatrick model (144)) are not yet entirely understood. This analogy with magnetic 
spin glasses has been very fruitful and in particular Hopfield made a breakthrough in the field of formal Neural 
Network by developing the analogy between a Neural Network with binary neurons and a spin glass. He showed that 
information can be stored in the minima of V and he proposed a method, inspired from the Hebb's rule to construct 
the interactions in a way such that the patterns to be learned correspond to the minima of V. 
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Indeed, assume that the J^j's are now given by equation (83), where the ^'^'s arc some patterns that we want to 
store and retrieve from the dynamics (87). For simplicity we take /Jfe = 1 in (83). We assume that the number of 
patterns , p, is < N. Since, in the limit g ^ go the phase space becomes {—1, 1}^, the takes binary values { — 1, 1}. 
Assume moreover that the vectors are mutually orthogonal, i.e. (C*^)^') = N5{k,i)- The Lyapunov function writes 
then: 



N p 
fe=l 

Note that in the case g ^ oo = sgn{ui). If x = then V = ^^^=1 i^lf i^jf ' TT E5=i EL2 ^UH]^' = 
^ Tf X]fc=2 i^^j^'^)'^ ■ Since the ^'^'s arc orthogonal one gets finally V = —N which is the absolute minimum 
of the magnetic energy. This results holds obviously for all patterns. Consequently, all patterns are absolute minima 
of V and they are stable^^. When the number of patterns is larger than A'', the patterns can not be all orthogonal. 
Therefore, the second term J2k=i (C^'^*^) pl^ys an important role, since it generates spurious minima of V. The 
exact analysis of the Hopfield model for a finite and infinite numbers of patterns have been performed by Amit, 
Gutfreund and Sompolinsky (2), in the thermodynamic limit A/' — > 00. For this, they uses spin glasses techniques 
such as the replica methods (in the case p — > 00). Their results have been rigorously proved in (67). 

When g is finite the neural network (87) has still the capacity to store and retrieve patterns via Hebb's rule (92). 

Moreover, for random symmetric Jij's the minima of V can still be computed by using techniques coming also from 
the spin glasses literature. For example, using a method developed by Bray and Moore (28) for spin glasses, Marcus 
& Westrevelt have been able to compute the number of fixed points for (87), when the size A — s- 00 (thermodynamic 
limit), in the case where J is a (symmetric) random matrices with independent entries and such that E[Jij = 0] and 
Var[Jij = jj]- They found equations for the fixed points which are similar to the Thouless- Anderson-Palmer equations 
(150) giving the mean-field solutions in the Sherrington-Kirckpatrick model, and the Bray & Moore techniques gave 
a similar result: the number of fixed points increases exponentially with the system size. 



C. Cooperative systems. 



When the synaptic couplings are not symmetric, (as in biological systems) there is in general no Lyapunov 
function, and many kinds of dynamics are possible. However, for some systems called cooperative systems one has 
still convergence properties, without Lyapunov function, but relying on a specific property of the flow, that preserves 
some pseudo order on the phase space. The results presented here are due to Hirsch (87). 

A dynamical system 



^=F,(u) (95) 



is called cooperative if: 



|^(u)>0, Vi^j (96) 

and it is called competitive if: 

|g(u)<0, y^^J (97) 

This has the following interpretation. As discussed in the section III.B the Jacobian matrix element |^ (u) measures 
in some sense the "influence" of the neuron j on the neuron i, when the system is in the state u. More precisely, it 



Actually, we have still to define the "limit" of the dynamical equations (87) when g — > 00. The Hopfield model with binary states uses 
a discrete time sequential dynamics (see (90) for details) 
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characterizes, to the first order in a Taylor expansion, the modification induced on Ui when Uj has a small variation. 
In cooperative systems the corresponding interaction graph has therefore only positive edges, whatever the state of 
the neural network and consequently, only positive circuits. They have moreover the following property. Assume that 
the phase space is convex and define the partial order u<v<^Mi<Wi, i = 1 . . .N. A cooperative flow preserves 
this order. Thus u(0) < v(0) ^ u(t) < v(f), Vt > (this corresponds to the positive feedback discussed above). 
Note that if F is competitive a reversal of the time arrow leads to u(0) ^ v(0) u{t) ^ v{t) > 0. From these 
inequalities, Hirsch (87) was able to prove that for a two dimensional cooperative dynamical system, any bounded 
trajectory converges to a fixed point. In larger dimension, one needs moreover a technical condition on the Jacobian 
matrix of F : it must be irreducible. Then Hirsch proved that the oj limit set of almost every bounded trajectory is 
made of fixed points. One does not have the same property for competitive systems (145). 

Some extension of these results have been recently made (74; 142). Though these works where intended to obtain 
mathematical results about metastability in the context of genetic networks, they hold in a very general context, 
and, in particidar, in the context of neural networks. In 1981, R. Thomas made the following conjectures (149). 1) A 
positive feedback loop in the graph of interactions of a differential dynamical system is a necessary condition for the 
existence of several equilibria. 2) A negative loop is a necessary condition for a stable periodic behavior. J.L. Gouze 
proved these conjectures under the hypothesis that the sign of the Jacobian matrix elements do not depend on the 
state. Consequently, the graph is the same everywhere in the phase space. The proof of the conjecture 1 has been 
extended by Soule in (142). The main idea in the proof of conjecture 1 is that if the dynamical system has several 
fixed points then F has several zeroes and thus cannot be injective. Thus knowing sufficient conditions for F to be 
injective, their negation give necessary conditions for F to have several zeros. The injectivity is ensured by conditions 
on the determinant of the Jacobian matrix. The proof of the second conjecture uses the fact that if all semi circuits 
(closed path in the non oriented graph) of length 2 < p < N are nonnegative then the dynamical system is equivalent, 
up to change of sign of some variables, to a cooperative system. Then there is no attracting periodic trajectory. 

As a conclusion, let us remark is that the notion of negative circuits is related to a notion of frustration introduced 
in the context of spin glasses. A loop is frustrated if the magnetic energy cannot be minimized for all the edges of 
this loop. This implies that the magnetic energy of this loop cannot reach its absolute minima and that sevcrals spin 
configurations lead to the same local minimum. The loop is frustrated because there are always "unsatisfied" edges 
where the spins are not in a configuration allowing them to minimize their energy. Flipping one spin may satisfy them 
but then others link will become unsatisfied. The notion of positive circuits is similar to the notion of non frustrated 
loop. Actually, one can obtain convergence results for symmetrically signed networks {sgn{Jij = sgn{Jji)) provided 
that the corresponding graph is not frustrated (24; 146). The notion of frustration has therefore nicely been adapted 
here. 

Finally, note that feedback effects can generate very complex situations, even if the dynamics is convergent. For 
example, in the case of symmetric synapses where a Lyapunov function exists, the presence of feedback terms (and 
frustration) induces a multiplicity of stable fixed points. This effect is analogous to the multiplicity of solutions for 
the Thouless- Anderson-Palmer (150) equations giving the various phases in the Sherrington-Kirckpatrick spin glass. 

D. Neural oscillators. 

What happens now if the synapses have no particular symmetry ? Actually, there are many possibilities including 
chaos. An example is presented below. But to end up the section V we would like to point out an easy way to generate 
oscillations in a system of two neurons having the dynamics (87). Consider therefore the dynamical system (87) with 
two neurons where, for simplicity, 9i = 02 = 0: 



where g controls the non linearity. It is easy to see that u = is always a fixed point. Moreover the Jacobian matrix 
at u = is: 



where / is the identity matrix and J the matrix of synapses. Therefore, the eigenvalues of DFq arc = — 1 + gsk 
where si, S2 are the eigenvalues of J^. We note therefore that the stability of the origin is determined, in this case, by 
the eigenvalues of J. We shall return back to this point in the section VI. The eigenvalues of J are straightforward to 



ill = —Ui + Jii tanh(gui) + J12 tanh(gu2) 
U2 = — W2 + ^21 tanh(gui) + J22 tanh(gu2) 



(98) 




(99) 
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compute and the eigenvalues of DFq are given by Xi^2 = — 1 — g '^^^\'''^'^ ± f — J22Y + 4Ji2 ^21- Consequently, it 

is possible to have oscillations in the system, provided that (Jn — J22)^ +4Ji2 J21 < 0. This imposes that J12 J21 < 0. 
Namely one neuron (say the first one) excites the second one while the second neuron inhibits the first one (see Fig. 
25). Note however that this is only a necessary condition. Actually, a sufficient condition corresponds to having a 
Hopf bifurcation destabilizing u = and generating stable oscillations. This happens whenever the two following 
conditions are fulfilled. 



(Jll-^22)'+4Ji2J21 <0 



(100) 
(101) 



For example, the following matrix (corresponding to the diagram drawn in Fig. 25) generates oscillations whenever 



-1 1 /■ 
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FIG. 25 Example of a system of two oscillating neurons (provided that g > 1). 



This example shows that once the synaptic symmetry is broken the dynamics usually does not settle onto a fixed 
point but on a more complex attractors. Even with two neurons the competition of excitation/inhibition can generate 
periodic oscillations. When the number of neurons increases, one can have pools of synchronized and oscillating 
neurons, as discussed in the section IV. F. But one can also have more complex situations ranging from oscillations 
to chaos. An example is given in the next section. 



VI. A COMPLETE EXAMPLE. 

In the previous section we have discussed various aspects of a recurrent model with convergent dynamics. Actually, 
in the early eighties a large effort was devoted to the study of convergent recurrent neural networks. Indeed, conver- 
gence was interpreted as a retrieval of a stored pattern (90). As shown above the symmetry of the synaptic connexions 
ensures, in the Cohen-Grossberg, model, the existence of a Lyapunov function and, consequently, the model is conver- 
gent. When the synapses are not symmetric, the dynamics can be quite a bit more complex, exhibiting a wide variety 
of dynamical regimes such as periodicity, quasi periodicity, chaos, existence of several complex attractors, etc... 

In the end of the eighties, some attention have been paid to recurrent neural networks exhibiting such complex 
dynamics. Indeed, the real brain is clearly a highly dynamical system and the convergence of the EEG to a fixed point 
is not a sign of a good health. From the biological point of view, accurate models have been designed to modelize 
temporal phenomena in the brain: synchronization of oscillations for feature linking (75), (59), (76), transition between 
coherent states (103), and chaos (15), (16), (17), (18), (68). Relying on recent neurophysiological results, the study of 
chaos in neural networks seemed a very promising way in at least two directions: the comprehension of the cognitive 
processes in the brain (139), (23), (42) and the development of new technologies involving the control of chaos and 
the massively parallel computability of neural networks. However, due to their complexity, they were very difficult to 
treat onto a mathematical ground and they lacked a theoretical background explaining the behavior of the networks 
in function of a few relevant control parameters. 

In this setting particularly astonishing experiments were made by Freeman and his collaborators on the olfactory 
bulb of the rabbit (63), (64). They suggested that the spontaneous dynamics of the olfactory bulb could be chaotic. 
But they also lead the authors to conclude that the recognition of a previously learned smell is manifested by a 
temporary reduction of the chaotic activity. On the basis of these experiments, Skarda and Freeman (139) proposed 
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an interpretation and a modeling scheme of the learning and recognition processes. In this scheme, the spontaneous 
dynamics of the neurons is chaotic and the retrieval of some previously learned pattern corresponds to a reduction of 
the chaotic attractor towards an attractor of lower dimension. During the alert waiting state, the network explores 
a large region of its phase space through chaotic dynamics. When the learned stimulus is presented, the dynamics is 
reduced and the system follows the lower dimensional attractor which has been created during the learning process. 

This idea is exciting but quite controversial, since, in particular, it is extremely difficult to measure quantities, such 
as fractal dimensions, on the basis of time series which are, intrinsically, non stationary. Nevertheless, this paradigm 
merits to be explored. For this, and to escape from the inherent limitations of data measurements in biological 
experiments, one possibility is to propose a model, that can be an oversimplification of the biological reality, but 
which captures some important features. The advantage of such a model is that one can simulate it numerically and 
have a better control on the time series. Also, sometimes, it is possible to obtain analytical results. 

This section is entirely devoted to such a model. Its structure is directly inspired by the Cohen-Grossberg model (but 
with a discrete time). Despite its simplicity, it displays an astonishing variety of dynamical behaviors and it has quite 
unexpected properties. Moreover, a relatively deep mathematical analysis can be performed combining concepts and 
methods from dynamical systems theory, statistical physics and ergodic theory. The subsections VI. A, VLB, VI. C, VI. D 
are devoted to a preliminary analysis of the spontaneous dynamics of the model, namely without learning. In the 
subsection VI. E we discuss its behavior when an Hebbian learning is applied and we show that a behavior similar to 
Freeman's paradigm is exhibited. Namely, this Neural Network has the ability to store information and to retrieve it 
by reduction of chaos. Finally, the subsection VI. F explores an important aspect and propose a new analytical tool 
to partially answer a basic question, arising naturally from the discussion above. Assuming that the model presented 
here as some relevance for brain dynamics, how can a chaotic network process information ? 

A. Model description. 

Let us consider a discrete time version of the equations (87): 

Xi{t+1) = f{Ui{t + l)) 

(102) 

Ui{t + 1) = J^3X3 (t) +0i i=l.-N 

where / is a sigmoidal function such as f{x) — tanh{gx) or f{x) — i±£££M2£i_ Henceforth, / maps IR to an interval 
[a,h] and the dynamics(102) occurs in a compact space f2 = [a,b]^. The parameter g controls the non linearity of 
/. This nonlinearity plays an important role in many aspects. Firstly, it is directly related to the transition to chaos 
described in section VI.C. But it has also an important influence when discussing the amplification/saturation effects 
on the propagation of a signal transiting via a neuron (see section VI.F). 

In this model, each neuron interact with each other (fully connected model). The "output" state Xi{t + 1) is a 
function of the weighted sum of the signals arriving at i at time t, Ui{t) = JijXj{t). We call Ui{t) the local field. 
Moreover, the "synapses" Jij are independent, identically distributed random variables, with expectation E[Jij] and 
variance Var[Jij] given by: 

mj] = ^ ; Var[Ji,] = ^ (103) 

such that the expectation and the variance of the "synaptic potential" JijXj{t) remains bounded as N ^ oo. For 
technical reasons we shall furthermore assume that the probability density of the Jjj's, p, obeys: 

(i) 3/3 > 1 s.t. / pl^{x)dx < oo 

(u) 3(5 > s.t. E [I Jf+'^l] < oo (104) 
(m) 3a > s.t. E [\J^-\] < n«",Vn > 2. 

Note that these conditions hold for a Gaussian or a uniform distribution. In the sequel the matrix of synaptic couplings 
will be denoted by J . 

In the section VI.E we shall discuss the effect of an Hebbian learning on the synapses and on the dynamical 
evolution. But at the present stage, assuming independence between the Jij 's may be viewed as an attempt to model 
a neural network initially "empty" of any information encoded in its synaptic structure [tabula rasa). Note that the 
synapses are therefore (almost-surely) asymmetric in this model. The situation is thus different from the previous 
section, where the symmetry allowed us to exhibit a Lyapunov function, ensuring convergence to fixed points. In 
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the present model, the attractors of the dynamics are in general not fixed points, but complex objects (e.g. strange 
attractors). 

In eq. (102) the quantities 9i play two different roles. In the present case (without learning) they correspond to 
a threshold in the neuron response. To take into account neuron diversity we assume that the Oi's are Gaussian, 
independent, and identically distributed random variables, such that: 

E[9i] = 9; VarlOi] = al (105) 

We call the vector {9i}f^i- 

In section VI.E we shall consider the effect of an input on the dynamics. In this cases, the 6i's will corre- 
spond to the input submitted to the neuron i. Also, we shall discuss in section VI.F the case where the input is 
time dependent. Finally, note that in eq. (102) the synapses and thresholds do not evolve in time (quenched disorder). 

The dynamical system (102) depends a priori on 'il* iV2 + AT + 1 parameters (AT^ synapses, N thresholds, and 
g). We call the vector A = (g, J, 0) the vector of microscopic parameters. A has therefore A^^ + N random entries. In 
the sequel, it will be useful to write the dynamical system (102) in the form: 

u(i + 1) = F [u(i); A] (106) 

Note that the Jacobian matrix as a simple form: 

D¥ [u; A] = JA(u) (107) 

where A(u) is the diagonal matrix such that A(u)ij = f'{uj)Sij. 

Let us now denote by Pj^g the joint probability distribution for couplings and thresholds, in a A' dimensional 

system. This probability distribution determines therefore the actual realization of the Jij's and ^j's. Moreover, Pj^g 

is determined by the parameters J, J, 6, erg. Hence, these parameters fix the statistical properties of the microscopic 
parameters and 9i. we shall call the parameters'^ g, J, J, 9, ag the macroscopic parameters. Note that we have only 



four independent macroscopic parameters because the dynamical system (102) is invariant under the transformation. 

g^gj. J^^^i^- Oi^j; (108) 



Hence J is somewhat irrelevant. 

From the dynamical system point of view widely developed in the previous sections, it is natural to seek the 
generic (in a topological sense) behavior of (102) when the microscopic parameters are varied. However, this is a 
formidable task and one may argue that, since these parameters are random, it is certainly more useful to investigate 
the generic behavior (in a probabilistic sense) when the macroscopic parameters are varying. In some sense, we 
substitute the analysis of the dynamical system (102), with uncountably infinitely many possible realizations of the 
microscopic parameters, by an "averaged" dynamical system depending on four independent deterministic macroscopic 
parameters. In this spirit, a few results are given in the next section, obtained by combining dynamical systems 
theory and probabilistic results about random matrices (section VLB). But, essentially, this approach is the core of 
the dynamic mean field theory, that will be fully developed in the chapter II. In the present chapter (section VI.D), 
we derive the mean field equations by an heuristic argument, and discuss their dynamical properties in relation with 
the dynamical system (102). 

Before entering into the detailed analysis let us make a last remark. The dynamics of the uncoupled neurons in the 
dynamical system (102) is rather poor. It indeed writes Xi{t + 1) = f{gJiiXi{t) + 9i). This dynamical system exhibit 
either a stable fixed point or bistability (appearing by a saddle node bifurcation) . Contrarily to the examples studied 
in the section II there is no Hopf bifurcation, no homoclinic loops, etc .... Nevertheless, the coupled system, as we 
shall see, has a rather rich dynamics. This provides a prominent example of emergent collective behavior. 



In this particular case g is both a microscopic and a macroscopic parameter. This is simply because all neurons have the same g. 
One can imagine a generalized version where the nonlinearity of the neuron i, gi, depends on the neuron and where the gi's are 
randomly distributed. In this case the gi would be additional microscopic parameters, while the parameters controlling their probability 
distribution would be additional macroscopic parameters. 
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B. Preliminary results 

Let us first establish a few preliminary results. Firstly, it is easy to show that, for each realization J{uj) of J , there 
exists a g value, gas{^), independent of 0, and given by: 



5as(c^) = ... (109) 

a\\J{uj)\\ 

such that F is a contraction whenever g < gas (33). In (109) || || is the operator norm induced by the euclidean norm 
and a is such that a = /'(O)- This result is straightforward since 

||F(u;A) -F(v;A)|| < sup ||i:'F(w; A)|| ||u - v|| = || ||u - v|| sup ||A(w)|| 
we SI weo 

where the last inequality holds from eq. (107). Since A(u) is a diagonal matrix such that Aij(u) = f'{uj)Sij and since 
/ is a sigmoidal function where the maximal slope is equal to ag, F is a contraction provided that a.g||i7|| < 1. The 
result follows. 

When F is a contraction the dynamical system (102) is absolutely stable i.e. it admits a unique fixed point, 
attracting all trajectories. The matrix is random and the result (109) holds for each realisation. One can obtain 
from this a statistical result by using a theorem proved by Geman (70). Provided that the J^ 's obey the condition 
(104(iii)), J' converges almost surely, when oo , to a finite value that can be explicitly computed (33; 35), 

depending on the parameters J, J. From this, one obtains the asymptotic limit of gas when N — > oo. This pro- 
vides, for finite A'^, an estimate of the g parameter values where the system is absolutely stable with a high probability. 



When g increases, one expects bifurcations leading to dynamical changes. When f{x) — tanh{gx) and when there 
are no thresholds, the function F(u; A) has the symmetry F(u; A) = — F(— u; A). Thus, u = is always a fixed point. 
Also, in this case L'F(O) — gj. Consequently, the stability of this fixed point is determined by the spectrum of 
the random matrix gj^ . Obviously, the eigenvalues of gj^ are proportional to the eigenvalues of with a coefficient 
g. J being real, the eigenvalues are either real or complex conjugated. When g increases, the spectrum is dilated 
and for sufficiently large g some eigenvalues are crossing the stability circle {z € C||2;| = 1} (see Fig. 26). However, 
the probability that several eigenvalues cross simultaneously this circle is zero^* if one excepts the case of a pair of 
complex conjugated eigenvalues. We expect therefore a destabilization of by a codimension one bifurcation. The 
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FIG. 26 Three possibilities of codimension one bifurcation occurring in (102) when tlie non linearity g increases. 



Having several eigenvalues crossing simultaneously the stability circle corresponds to impose algebraic relations of codimension larger 
than one between the coefficients of the characteristic polynomial of Since the Jij's are selected randomly the probability to fulfill 
these algebraic conditions is zero. 
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possible codimension one bifurcations for a map having the symmetry F(u; A) = — F(— u; A) are described in Fig. 26. 
They are (see the appendix for details): 

• Case I. Pitchfork bifurcation. The fixed point destabilizes and two symmetric stable fixed points appear. 

• Case II. Period doubling bifurcation. The fixed point destabilizes and stable periodic orbit of period two appears. 

• Case III. (Discrete time) Hopf bifurcation. The fixed point destabilizes and stable periodic orbit appears. Note 
that, stricto-sensu, orbits of period 2,3 and 4 do not correspond to a Hopf bifurcation (the normal form is 
different, see (13) for details). Orbits of period 3 and 4 are observed for small iV's (57). 

Call /o(J') the spectral radius of J (the value of the largest modulus of the eigenvalues). Then the destabilization 
occurs when: 



(110) 

This is a random variable. However, the statistical behavior of random matrices obeying the conditions (104 (i), 
(ii)) is well known when the size tends to infinity (73). The limiting spectral density converges almost surely to the 
uniform density in the disc of center and radius J in C. Consequently, go converges almost surely to j and the 
destabilization value is given by go>^ = 1- Note that the same result can be obtained from the dynamical mean field 
theory (see (148) and (47) for a continuous time version of (102)). 

The repartition of eigenvalues is also known in the finite size case (60). One can show that there is an over density 
of real eigenvalues that disappear in the limit N ^ oo. Consequently, for finite size, one observes destabilization 
by pitchfork and flip bifurcations, but the Hopf bifurcation becomes more and more frequent when N increases 
(32; 36; 57). Finally, in the infinite system an infinite number of eigenvalues cross simultaneously the unit circle. This 
corresponds to a sharp transition from fixed point to white noise discussed in section VI.D (see also (36)). 

Let us now make a remark about the Hopf bifurcation. As we somehow anticipated in the sections IV.F,V.D 
oscillations arise because there is a competition between excitation/inhibition effects among the neurons. Actually, 
one expects from the study performed in IV. F to have, near the bifurcation, pools of almost synchronized neurons 
oscillating coherently. This is revealed in the study of the correlation function which has usually a bloc structure as 
revealed from example in (51). Note also that the period of oscillations is generically irrational. Finally, the results 
above the spectrum of imply that the phase v of the largest eigenvalue, generating the Hopf bifurcation, is uniformly 
distributed between [0+,7r]. Hence Pro6[0+ < u < 6] = ^ and, since the period is T = ^, Prob[T > r = ^] = |. 
Therefore the probability density of the period is Pt(t) = Thus, there is a high probability to have oscillations 
with a low period. 

These results have been obtained by combining elementary results from dynamical system theory, holding for each 
realization of the disorder, and convergence results in random matrices theory. The convergence results, holding 
when N oo, are then used as a guideline for a typical realization of the finite dynamical system. They are 
however quite restricted. For example, we have assumed that the system has the symmetry F(u; A) = — F(— u;A). 
But when we consider the equation (102) with thresholds, this symmetry disappears. Then, the fixed point of the 
absolutely stable regime is a random variable. Moreover, when g > gas new fixed points can appear by saddle-node 
bifurcations: they are also random. Finally, we have been able to analyze the first bifurcation relatively easily but, 
after the destabilization the usual techniques (central manifold reduction, normal forms) are difficult to handle since 
the coefficients are random (hence, for example, the eigenvectors of are random). 

One has therefore to develop an alternative statistical approach. This is done in the section VI.D. Before this, we 
discuss in the next section the typical behavior of the dynamical system (102) when g further increases. The results 
presented are a combination of genericity results in dynamical systems theory and numerical simulations. 

C. Transition to chaos 

Numerical simulation is a fundamental tool for the exploration of the wide dynamical richness of the model (102). 
But clearly, exploring the parameters space of this system at "random" , without any preliminary idea of what is going 
on is like "searching a needle in a straw pile" (it is in fact a bit more tricky since a straw is only a three dimensional 
object) . Indeed, basically, dynamical systems are structurally stable on wide ranges of parameter values and only the 
points where structural stability fails (bifurcations points) matter. But bifurcations occur for parameter values usually 
located on manifold of smaller dimension than the ambient space. For example the codimension one bifurcations 
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discussed above correspond to a, J\f — 1 — N'^ + N manifold in the microscopic parameters space. Since we select the 
Jij's and 9i^s with an absolutely continuous probability distribution (i.e. having a density), the probability to fall on 
a bifurcation point is zero. Obviously, since we are seeking statistical properties, we are rather interested in locating 
the bifurcation points in the macroscopic parameter space. Having a bifurcation map in this space would correspond 
to having statements such as: "If you fix the parameters J, 9, ag in this region of the macroscopic parameters space, 
and if you vary g between such and such value, then, typically you will observe this type of bifurcation" . 

For this, we need to have theoretical guidelines. The preliminary results given above are an example of such 
guidelines. The mean field approach briefly discussed below provides additional hints. Consequently, the numerical 
simulations described in the present section have been made with the informations given by these theoretical results, 
plus a few standard and generic facts in bifurcation theory. These facts are: 

• Breaking the symmetry F(u; A) = — F(— u; A) transforms pitchfork bifurcations into saddle-node bifurcations as 
depicted in Fig. 27. We observe indeed such bifurcations and we have an analytical way to locate them (see the 
next section). 




FIG. 27 Effect of breaking the symmetry F(u; A) = — F(— u; A) on a pitchfork bifurcation. 

• The fixed points can be destabilized when g increases. They generically do it by Hopf bifurcation (namely with 
an increasing probability as N increases). 

• As shown in Fig. 28, after the first Hopf bifurcation, the standard scenario is the "Ruelle-Takens-Newhouse" 
transition to chaos by quasi periodicity (130), (though our system is a discrete time system). As g increases 
the limit cycle generated from the first Hopf bifurcation destabilizes by a second Hopf bifurcation giving rise to 
a two dimensional (T2) torus. Near the bifurcation, the trajectories densely fill the torus since the frequencies 
corresponding to the first and second Hopf bifurcation are, in general, irrational. However, a further increase of g 
leads to a frequency locking: the frequencies corresponding to the first and second Hopf bifurcation synchronize in 
a rational fashion and the trajectories are periodic orbits on the torus. Though frequency locking is structurally 
stable, increasing enough g finally lead to chaos, by different ways (for a detailed explanation in general models 
see (69; IIO); for a detailed description of this model see (35)). Note however that there may exist "re- 
stabilisation phases" when g further increases. This corresponds usually to the crossing of "Arnold tongues" 
where the dynamics locks on a quasi periodic orbit. An example is given in Fig. 30a where have plotted the 
first and the second Lyapunov exponents. The first Lyapunov exponent increases with g except at some points 
where it takes a zero value. Since the second Lyapunov exponent is also zero this corresponding to a reduction 
of the chaotic dynamics on a T2 torus. If one excepts these points, the positive Lyapunov exponents and the 
fractal dimension of the strange attractor increases as g increases. 

In fig. 30b we have plotted the Lyapunov spectrum for g = 3.5. One notes that, in the example chosen, there 
is only one positive exponent. Thus the corresponding (Kaplan- Yorke) dimension is low {Dky — 1.967). More 
generally, one observes that the strange attractor is usually a low dimensional object (compared to the dimension 
of the embedding space). One consequence is that an arbitrary perturbation of a point on the attractor has 
almost all its components outside the attractor. Note finally that this transition to chaos generates resonances 
peaks in the power spectrum (Fig. 29) some of them resulting from the Hopf bifurcations. Thus, even if in the 
chaotic regime the power spectrum is continuous, it is not flat, like white noise, but it has peaks or resonances. 
These remarks lead to important issues discussed in the section VI. F. 

• As A'^ increases the transition to chaos occurs on a g range becoming more and more narrow. This leads to the 
conjecture that a sharp transition from fixed point to infinite dimensional chaos occurs in the thermodynamic 
limit. This conjecture is related to the observation above that the eigenvalues of the Jacobian matrix accumulate 
on the stability circle as A'' ^ oo. Exact equations and analytical description are discussed in the section VI. D. 
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FIG. 28 An example of transition to chaos by quasi periodicity in the model (102). We used the representation m{t + 1) versus Tn{t) 
where m{t) = -^X^iLi^il*) is the empirical average of the neurons state at time t. This representation provides a projection of the 
dynamics of m{t) in a two dimensional phase space. The insets represent the evolution of m{t) versus t. In the two last figure Ai is the 
largest Lyapunov exponent. 



D. The mean-field dynamical system. 

The mean field "approximation" is quite well known in statistical physics. Though it is stricto-sensu wrong in many 
cases (for example it gives a wrong critical temperature for the Ising model), it provides often an astonishingly good 
qualitative insight in the description of many models of phase transitions in statistical physics. Also, for some models 
(such as the Curie- Weiss model) it is exact in the thermodynamic limit. In the field of neural networks the use of 
mean field approaches has a long history, for analogic networks (148) but also for spiking networks (see for example 
(29)). The chapter II is entirely devoted to mean field approaches, and in the present chapter, we shall focus only on 
the model (102). 

Basically, the mean field approach applied to this model consists in assuming that the a;i(t)'s are independent 
from each others and independent from the J^'s !! Though this looks very rough, this approach leads to exact 
results that can be rigorously proved (see chapter II). It can also be justified at an heuristic level (36): one easily 
shows that the key ingredients ensuring the success of this approach are the independence of the Jij 's and the fully 
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FIG. 29 Power spectrum of the neuron in the transition to chaos corresponding to the previous figure. 



connected structure of the model. Hence, the mean field approach breaks down as soon as some correlation between 
the Jjj's exist (e.g. after learning). Note also that it breaks down when the J^'s are symmetric. More precisely, one 
has to correct the mean field equations derived below by adding a feedback term corresponding to the delayed ac- 
tion that a neuron has on itself (36; 46; 114; 150) (this action cancels, on average, when the 's are independent (36)). 

Assume therefore that Xi{t)'s are independent from each others and independent from the Jij's. Then the central 
limit theorem states that the "local fields" Ui{t + 1) = J2j Jij^ji^) + Gaussian processes in the limit N ^ oo. 

Moreover, they are independent and identically distributed. Hence the joint probability of the Ui(i)'s factorizes in 
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FIG. 30 30a Evolution of the first two Lyapunov exponents when g increases. 30b Lyapunov spectrum in the chaotic regime when g = 3.5. 



an (infinite) product of identical Gaussian distributions. To cliaracterize tire Gaussian distribution at time t one 
needs the average value /^(t) = E[ui{t)], the variance v{t) = E[u'j{t)] — and the time covariance A{t,t') — 

E[ui{t)ui{t')] — ix{t)ix{t') (note that the left hand side terms are independent of i since all the local fields have the 
same distribution). It is straightforward to see that these quantities are functions of m{t) = E[x{t)], q{t) = E[x^{t), 
and C{t,t') = E[x{t)x{t')] (see eq. (Ill), (112), (113)) below. Finally, since x{t) = f{u{t)) and since u{t) is Gaussian 
one obtains: 



m{t) 



Jm{t) 

+ 00 



/2tt 



(111) 



v{t + l) 



rq{t) + ae 

+ 00 



'2ti 



(112) 



f{hy^ + fl{t))dh 



A{t + l,t' + 1) = J^C{t,t') + a'^g 



Clt.t') = 
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(113) 

Kt) ) / (/I'T^ + n{t')) dhdh' 



These are the dynamic mean field equations of the model (102). The parameter m,q are called "order parameters" 
in the statistical physics literature. They characterize the emergent behavior of a system with a large number of 
degree of freedom and they exhibit drastic changes corresponding, in statistical physics, to phase transitions, and in 
our context to a macroscopic bifurcation. 

Let us now make a few remarks. Firstly, these equations can also be derived from the computation of a generating 
functional for the probability distribution of the trajectories {u{t)}^^. This has been developed by Crisanti et al. 
(47) for a continuous time version (without thresholds) and by Molgedey et al. (116) for a discrete time version. 
Both approaches lead obviously to the same equations. But they also deal with the same type of convergence 
namely weak convergence. As said above, the idea below the mean field approach is to have informations about 
the "average" behavior of the dynamical system (102). This is what we have obtained, but in a very rough 
sense. The equations (111), (112), (113) tell us about the evolution of the average value of u{t) when the average is 
performed over infinitely many realizations of the disorder. But, weak convergence does not give any information 
about one typical system whose size tends to infinity. For this, one needs a stronger convergence, the almost- 
sure convergence^^. The large deviations approach developed in chapter II will allow us to show the almost sure 



Almost sure convergence corresponds to the statistical physics notion of self averaging. The empirical average of a quantity in one 
realization of the disorder converges with probability one to the average of this quantity over the disorder 
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convergence (and incidentally that the equations (111), (112), (113), derived from a "questionable" Ansatz, are correct). 

Let us now discuss these equations, their solutions and their interpretation. One remarks firstly that, for t = t' , the 
equation of A{t, t) is the equation of (112) for the variance v{t) (as expected since A(f, t) = v{t)). Also, A(f, t') < v{t). 
One can therefore write the equation (113) in the form 

A(t + 1, i' + 1) = i/,, j,^,,,.. [A(t, t')] j'C{t, t') + (114) 
where H is defined only when A^(t, t') < v(t}v{t') and is given by (113). 



Having these equations in hand, the idea is now to study the reduced dynamical system (111), (112), (113) and to 
infer informations about the typical dynamics of (102). More precisely, we are interested in the time asymptotic that 
corresponds to a stationary regime of (111), (112), (113). The stationary equations are given by: 

l_i = Jm + e (115) 

/•+°°e-^ / - 

m = —j=f{hJj'^q + al + Jm + e)dh 

J — oo V Stt 



V = rq 

r'+oo 



(116) 



/ —j=f{hJj'^q + ul + Jm + e)dh 



A(t - f ) = fiC(t - t') + a» = ff,,.7„,,s,„ |A(« - «')] (117) 

J_oo J-oo V27r V2^ \ J 

These equations give important informations about the statistical behavior of the model (102) with an increasing 
accuracy when the size increases. Moreover, m, q, A act as "order parameters" allowing us to distinguish different 
dynamical regimes and to draw an effective bifurcation map in the space of the macroscopic parameters. Let us list 
a few results (32), (34), (35), (36). 

1. In the absolutely stable regime the fixed point is a random variable. One checks numerically that the corre- 
sponding distribution of the local fields is Gaussian and the equations (115), (116) give the mean and variance 
of this distribution with a very good accuracy. For fixed N the empirical mean and variance computed over a 

large number of networks is close to the theoretical values. Moreover, the statistical dispersion of these empirical 
values decreases as N grows (as expected from the almost-sure convergence proved in Chapter II) . 

2. The equations (115), (116), have several solutions in some regions of the macroscopic parameter space. More 
precisely, they exhibit saddle-node bifurcations. The critical values where saddle-node bifurcations occur in the 

space of macroscopic parameters can be computed from equations (115), (116). It is remarkable that these 
bifurcations have a direct correspondence with the saddle-node bifurcations observed in the system (102) in the 
following sense. 

On one hand if one fixes the parameters ^,cre, J in a region where the mean field equations predict a saddle- 
node bifurcation as g increases one observes indeed (in general) saddle-node bifurcations in the system (102). 
Of course, the exact g value where the bifurcation occurs is random and depends of the actual realization of the 

disorder. However, if one performs a statistical analysis of these values, one finds an average value close to the 
value predicted by the mean field equations. Moreover, the empirical variance decreases with the system size. 

On the other hand, the various fixed points appearing from saddle-node bifurcations in the dynamical system 
(102) are also random. But the coordinates of these points (in the Ui space) are distributed according to a 
Gaussian distribution whose mean and variance are in good agreement with the value obtained from the fixed 
points of (115), (116). 

As a conclusion, the analysis of the fixed points of (115), (116) allow us to draw a bifurcation map in the 

macroscopic parameter space giving the average g value where saddle-node bifurcations occur, for a given value 
of the parameters 9,a0,J. It gives also the probability distribution of the corresponding fixed point in the 
dynamical system (102). 
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3. Once we know the statistical distribution of the fixed points, one can compute a destabilization condition by 
estimating the spectral radius of the Jacobian matrix in the same way as we did above (but now the distribution 
of eigenvalues depends on the distribution of the fixed points (34)). This condition is given by: 



5op(JA(x*)) = 1 (118) 

where p{) is the spectral radius, x* the fixed point, and A the diagonal matrix Aij(x*) = f'{u*)Sij. This 
condition generalizes (110) since for f{u) = tanh{gu), x* = and A(0) is the identity matrix. From this one 
obtains the average g value where destabilization occurs, for a given value of the parameters 9, ao, J. Moreover, 
the Jacobian matrix as similar spectral properties as in the case F(u: A) = — F(— u; A) and when a fixed point 
destabilizes it does this (generically) by a Hopf bifurcation. Therefore the mean field equations allow us to draw 
a bifurcation map in the macroscopic parameter space giving the average g value where a Hopf bifurcations 
occurs (see ref. (32)). 

4. A finer analysis of the complete set of equations ((115), (116), (117)), and especially of the equation for the time 
covariance (114)), reveals that there are in fact two regimes. The equation (117) admits always the solution 
A = V. This solution is stable for the map (114) if ^ < 1. The destabilization condition is therefore given 
by20: 

^(^A = v) = J^J ^-j=r{h^J^q + <jl + Jm + e)dh = l (119) 

This equation defines in the space of macroscopic parameters (g,J,9,a$) a codimension 1 manifold which 
separates this space into two regions. 

In the region corresponding to ^ < 1 the solution A = w is stable and it is the only solution of (117). The 
asymptotic stochastic process described by the mean field equations is then stochastically equivalent to the 
Gaussian process given by: 

y{o) = X 

Y{t + l) = Y{t) 

where X is a Gaussian random variable M{ijl,v). Henceforth, F is a process with almost-surely constant 
trajectories. Its interpretation is easy: it corresponds to a regime of (102) where we have only fixed points. 

In the other region, one can write v5{t,t') + A* {1 — 5{t,t')) = {v — A*)6{t,t') + A* . Consequently, the asymptotic 
stochastic process described by the mean field equation is stochastically equivalent to the Gaussian process given 
by: 

^(0) = X 
Z{t+1) = Z{t)+B{t) ^^^^^ 

where X is a Gaussian random variable .A/'(/i, A*) and where B{t) is a white noise with zero mean and variance 
{v — A*). Z{t) is therefore the superimposition of a process with almost-surely constant trajectories plus a white 
noise. 

It is also remarkable that the equation (119) corresponds exactly to the equation of destabilization of the fixed 
point. We conclude therefore that the crossing of the manifold (119) corresponds, in the infinite system, to a 
sharp transition from fixed point to infinite dimensional chaos. Note however that this "manifold" is a very rough 
representation of the edge of chaos for finite size systems. Indeed, it is known (69; 110) that in the transition to 
chaos by quasi periodicity, the edge of chaos has a fractal structure corresponding to the intersections of Arnold 
tongues. 



The equations ((115), (116)) arc similar to the Slicrriiigtoii-Kirkpatrick equations describing the Sherringtou-Kirkpatrick spin glass 
model(144) at high temperature, while the equation ((53)) corresponds to the De Almeida-Thouless line. A detailed discussion of this 
aspect has been done in(34; 36). 
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FIG. 31 An example of bifurcation map. The surface drawn in the parameter space 6,<7g,g corresponds to the sharp transition from fixed 
point to chaos, obtained from the mean field equations in the thermodynamic limit. 



The theoretical results described in the sections give a fairly good description of the various dynamical regimes 
generically exhibited by (102). However, mean field equations have the drawback to hold only when the size of the 
system tends to infinity. And we have just seen that this limit is rather poor (either fixed points or white noise). 
Therefore, though mean field equations are a good guideline for describing the statistical behavior of (102) they miss 
a lot of important aspects: intermediate regimes between fixed points and chaos, dynamical properties of a given 
realization of the network, etc. In the next section we depart from the rough vision provided by the mean field theory 
and develop two aspects drastically related to the finite size system. 

E. Hebbian learning effects. 

Let us now consider the effect of Hebbian learning on the dynamical system (102). For this, we return back to the 
recipes discussed in the section III.C. Learning is based on the modification of synaptic connections between neurons. 
In the present context, this is interpreted as a a slow evolution of the synapses Jij when the network is submitted to 
a pattern that one would like to "teach" to the network. In our model, a pattern is a vector ^ = {^i}i=i..]v and a 
presentation consists in adding the vector ^ to the vector of thresholds 9 (i.e. Oi ^ 9i + i = 1 . . . N). 

Several (many) learning rules can be proposed, based on the recipes presented in the section III.C (see for example 
(51)). A straightforward implementation, very similar to equation (82) is: 

Jij{t + l) = XJij{t) + ^{xi{t + l)-r])x{xj{t)-r])xH{xj{t)-r]) i,j = l..N; T>1 (122) 

The parameter < A < 1 corresponds to a decay of the synapse when it is not used. H is the Heaviside function 

{H{x) = 1 if .T > and otherwise), rj defines a level of activity allowing us to decide whether a neuron is "active" 
at time t {x{t) > rj) or "silent" {x{t) < rj). Consequently, the term {xi{t + 1) — r]){xj{t) — ri)H{xj{t) — rj) corresponds 
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to modifying the synapse only if the post- synaptic neuron is active. This corresponds to the fourth^^ "recipe" in 
the section III.C. Finally, a weight cannot change its sign (this corresponds to demanding that a synapse cannot 
switch from excitatory to inhibitory or vice- versa) . 

On biological grounds, the learning rule (122) can be interpreted as follows. The synaptic weight Jij con- 
nects the neuron j to the neuron i and the output signal emitted by j at time t is transmitted to the neuron 
i at the next time step with the weight Jij. Let us assume that A = 1 (no forgetting). Then, the learning 
rule has the effect of enhancing the synaptic strength Jij if the neuron j is active at time t and if the neuron i 
is active at time t+1. On the other hand, if j is active at time t and the neuron i is inactive then the synapse decreases. 

The joint evolution of (102) and (122) occurs as follows. The initial couplings and thresholds J^^ , 9^ are fixed to an 
initial random value determined by the probability distributions (103), (105), for a determined value of the macroscopic 
parameters (J,^, Cg). These parameters and g are fixed such that the corresponding dynamics is chaotic. The values 
of these parameters can be roughly determined from the bifurcation map described above (see Fig. 31). 

After a sufficiently long time such that the neurons dynamics has "reached" its chaotic attractor, one presents a 
pattern ^. This means that one modifies the thresholds: 9} — 9^ ^i. The weight Jij are not modified at this stage. 
The pattern is a random vector ^ whose entries are independent, identically distributed, Gaussian, with a mean ^ 
and a variance cr^. Henceforth, each neuron feels an effective threshold 9j = 9i + ^i. This modifies the dynamics. 
However, from the macroscopic parameters point of view, this amounts to have the transformation 9 9 -\- ^, 

^ Hence, it is still possible to know the average behavior of the perturbed system by using the bifurcation 

map. In the experiments described below, the pattern is chosen such that the perturbed dynamics remains chaotic. 
Then one iterates the learning procedure (122). The stimulus is always present. Once the learning phase is finished 
one removes the stimulus ^ (i.e. the thresholds are reset to their initial value 0°). 

What is the effect of the learning rule (122) (1) on the neurons dynamics ? (2) on the synapses ? The typical effects 
on the neurons is depicted in Fig. 32. One observes generically an inverse quasi periodicity route. Namely, starting 
from a chaotic attractor, the modification of the Jij's by the Hebbian rule (122) leads first to a torus, then to 
a limit cycle and, finally, to a fixed point (with possibly a crossing of several Arnold tongues leading to temporary 
synchronizations). Thus, too long a learning phase basically "kills" the dynamical activity (see Fig. 32). 

■I.'., 1 . . 1 1 . . 1 




FIG. 32 Fig. 32a Inverse quasi periodicity route induced by learning. In this simulation A'^ = 64, g = 8 a = 10 = 0.5, A = 1 (no 

forgetting). 30000 learning steps are represented. The plotted quantity is the average value of the output states m{t) = jf ^i^^i ^iit). 
Fig. 32b Graphical representation of the learned pattern. 

Now, assume that we stop learning when the systems is in an intermediate phase (e.g. quasi periodic or periodic). 
Different results are possible depending on the time where we stop learning but also on the pattern, the spontaneous 



We have removed the condition that the pre- synaptic neuron is active. Indeed, in this naive model, a term {xi{t + 1) — r]) X H(xi{t + 
1) — r]) X {xj{t) — Tj) X Hixjit) — Tj), always positive or zero, would lead to an increase of the synapses linking active neurons and to an 
exponential decay of the other synapses toward 0. Hence we would rapidly have a network composed by positive synapses only, with a 
value increasing in time. Thus, rapidly, all active neurons would saturate. 
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dynamics, the learning rule etc .... Nevertheless, it is possible to observe the following phenomenon, reminiscent of 
PVeeman's paradigm. In some cases, removing the pattern when the activity of the network+pattern is periodic gives 
back a strange attractor. Then, a new presentation of the pattern leads back to the limit cycle. An example is given 
in Fig. 33. The initial regime is chaotic (Fig. 33a I) and presenting the pattern does not change the chaotic nature 
of the dynamics (Fig. 33a II). Obviously, this changes the attractor, but nothing significant happens. In particular 
a glance to Fig. 33b does not reveal any clear cut effect induced by the pattern presentation, before learning. The 
situation is drastically different after learning. If one stops the learning phase corresponding to the figure 32 after 
11.000 learning steps one ends with a periodic attractor (Fig. 33a III). Then, removing the pattern leads back to 
chaos (Fig. 33a IV). Again, the form of the attractor is different from (Fig. 33a I,II) but observing the dynamics does 
not tell us that learning has been performed. However, a new presentation of the pattern induces a sharp reduction of 
the dynamics onto the limit cycle (Fig. 33b). Since, the pattern presentation didn't lead to this effect before learning, 
this property has been acquired via learning. 

Consequently, in this situation, the learning process associates to a given pattern a dynamical pattern, and recog- 
nition is manifested by a dynamical reduction from chaos to the associated dynamical pattern. We have therefore 
another possible interpretation for the loose statement given in the section III.C: "a pattern is memorized if the 
neural networks has acquired, via learning, the capacity to dynamically evolve towards a "state" "associated to the 
pattern", provided that it was "suitably prepared" ". Here the "state" is an attractor^^, and "suitably prepared" 
means that we add the pattern ^ to the vector of thresholds 6 (pattern presentation). Hence, the effect of learning is 
quite different from the Hebb-Hopficld learning where a pattern is associated to a fixed point and "suitably prepared" 
means choosing an initial condition in the attraction basin of the pattern. 




FIG. 33 Learning and effect of a pattern presentation after learning. Fig. 33a Attractors. I. Attractor before learning and before the 
pattern presentation. II . Attractor before learning when the pattern is presented. III. Attractor after 11000 learning steps with the 
pattern. IV. Attractor after 3000 learning steps without the pattern. Fig. 33b Time trajectories. 

The remarkable fact is that the learning dynamics has lead the system in a state different from the initial one. 
Without excitation by the stimulus, the neuron dynamics is chaotic and there is no apparent difference between this 
case and the situation before learning. More precisely, certainly the learning phase has changed the characteristics of 
the strange attractor, but this change does not tell us anything about the fact that an information has been encoded 
in the network. This fact is revealed only if one presents the stimulus and its manifestation is drastic (remember that 
the presentation of the pattern before learning didn't reduce the dynamic) . 

This observation raises however many questions in particular with respect to the robustness of this behavior, and 
the mechanisms leading to it. We postpone these questions to the end of the section and we investigate now the 
second point listed above. What is the effect of the learning rule (122) on the synapses ? 

The remarkable fact is that no clear cut changes are observed even if the learning phase is long. Obviously the 
Jij's are modified by the learning rule (122) but there is no striking change in the structure of the matrix or in 



To be more precise the state is an ergodic probability measure with support on the attractor. A natural choice is the SRB measure 
introduced in the appendix. Thus, in the present context, the notion of state is closer to statistical mechanics framework where a 
(macro)state is a probability measure on the phase space. 
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the histogram of the Jij's, even if these infinitesimal changes in the J^j's are sufficient to modify the dynamics. An 
example is given in Fig. 34. After 11000 times steps, the dynamics settle onto a limit cycle but the histograms and 
the matrix looks very much like the initial one. To observe significant changes one has to iterate the learning phase 
far beyond the time where the dynamics has died. The Fig. 34 shows the distribution of the J^ 's and the matrix after 
10^ time steps. Here a clear modification is revealed. The weights emitted by some neurons have increased, while 
the others have not been modified. But the time scale to observe a significant change is substantially larger than the 
time necessary to have a dynamical reduction. 




FIG. 34 Fig. 34 Effect of learning on tfie synapses. 34a Histogram of the Jij's. I Initial. II After 11000 learning steps. III. After 10^ 
learning steps. 34b Matrix ^7 I Initial. II After 11000 learning steps. III. After 10^ learning steps. The radius of the circle is proportional 
to the absolute value of the synapses. Blue circles correspond to inhibitory synapses and red circles to excitatory synapses. 

From the theoretical point of view, one is far from the degree of understanding of the model without learning and 
there is no quantitative theory allowing to predict and control the effect of learning. As a matter of fact, the dynamic 
mean field cannot be applied, since the learning dynamics (122) creates correlations between the weights. However, 
one can give the following heuristic explanation of the phenomenon. First, the mean field approach developed in 
the section VLD has left us with a somewhat misleading picture of the neural network. Indeed, in the mean field 
treatment all neurons are equivalent and thus they have all the same level of activity. This is correct if one considers 
the activity of the neurons averaged over a large number of networks. But the situation is different when one considers 
one particular realisation of the Jij 's. In the figure 35 we have represented the time averaged (see the appendix) value 
of Ui and Xi, in the various phases of the learning procedure. In the first row we have represented the average output 
activity {xi) only for the neurons such that {xi) < r/ = ^. Thus these neurons are (on average) active neurons. 
Though the learning rule (122) uses the instantaneous activity of the neurons and not the average (for a variant of 
this rule, see eq. (123) below), this representation gives us an indication of the repartition of "active" and "silent" 
neurons. This repartition is clearly not uniform, since it results from the interplay of the neuronal connections Jij's 
and the non linearity of the transfer function (this interplay and the resulting properties are discussed in more details 
in the section VI. F). 

The pattern presentation as a direct but weak effect on the local fields and an even weaker effect on the activity 
(the pattern is represented in Fig. 122b; note that obviously (xi) — (/(uj)) /((""«)))• The learning rule selects 
then the active neurons and modify their outgoing synapses in the following way. Assume that j is an active neuron. 
Then if i is active Jy increases. Thus Jij becomes more excitatory if it positive and it becomes less inhibitory if it 
is negative. On the opposite, if i is "silent" then Jij decreases. Thus Jij becomes less excitatory if it positive and it 
becomes more inhibitory if it is negative. In all other cases Jij stays constant (for A = 1). If we admit that one step of 
learning has a small influence on the level of activity of the neurons^^ then the picture remains essentially the same at 
the next learning step. Thus, in this rough picture, we have a set of active neurons whose outgoing synapse gradually 
evolves. The excitatory links towards active neurons become more and more excitatory, the inhibitory links towards 
silent neurons become more and more inhibitory; in the same time the excitatory (inhibitory) links towards silent 



This can be assumed away from bifurcations point (see section IV) but it is incorrect near a bifurcation. 
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FIG. 35 Fig. 35 Effect of learning on the activity of the neurons. First row. Average output activity (xi) of the neurons Second row. 
Average value of the local field 



(active) neurons decay to zero and eventually vanish since a weight cannot change its sign. Consequently, a very 
long learning phase leads to an histogram such has Fig. 34a III, with a high peak at zero, two bumps corresponding 
to excitatory and inhibitory synapses with a large^** absolute value, and finally a background of synapses that have 
essentially not been modified during learning. The active neurons become "hubs" for the dynamics, in the sense that 
they have a relatively large connectivity and some weights with big absolute values. This corresponds to the vertical 
bands with big circles revealed in Fig. 34c. One also has horizontal bands with mainly either red or blue big circles. 
The "red lines" corresponds to the links received by active neurons coming from active neurons, while the blue ones 
corresponds to the links received by silent neurons also coming from active neurons. It is thus remarkable that the 
Hebbian like learning rule (122) leads to a structuration of the network '^^ into "pools" of neurons. Finally, from 
the dynamical point of view, since active neurons become more and saturated the dynamics converges to a fixed point. 

This pictures gives us a fair understanding of the (somewhat trivial) behavior of the system when learning is 
performed on long time scales. But, what about the small time scales and what about the inverse quasi periodicity 
route ? For this, let us use the wisdom acquired in the preceding subsections. The dynamical system (102) can be 
represented by a (randomly) chosen point in a space of parameters with J\f = N"^ + A'^ + 1 dimensions. In this space, 
many "critical" manifolds exist, whose crossing corresponds to various type of bifurcations. As discussed above a 
complete investigations of this space is impossible but standard results in dynamical systems theory, completed with 
numerical simulations and mean field theory have allowed us to roughly locate the "boundary of chaos" as a function 
of the macroscopic parameters. Note however that the "bifurcation manifold" obtained from the mean field approach 
in the figure 31 is a very rough representation of the edge of chaos. Indeed, it is known (69; 110) that in the transition 
to chaos by quasi periodicity, the edge of chaos has a fractal structure corresponding to the intersections of Arnold 
tongues. Thus the transition is usually not sharp when one modifies the parameters but one has succession of phase 
locking with various rotation numbers and chaos (see e.g. Fig. 30a and 32). 

On the other hand, the manifold corresponding to the destabilization of the fixed point has a nicer structure. It is 
indeed given by eq. (118) (7op(i7A(x*)) = 1. Now, the learning dynamics corresponds to a motion of the representative 
point of the dynamical system in the subspace of synaptic weights, while the presentation or removal of the pattern 
correspond to a translation in the subspace of thresholds. These motions lead to bifurcations when crossing critical 
manifolds. Consider now the destabilization condition (118). Since learning has the effect of slowly increasing the 



Note that there is no upper or lower bound on the synapses in the learning rule (for a variant see eq. (123)). Thus, the modified 
synapses diverge asymptotically. 

The Hebbian learning generates in fact small word structures, as shown in (21). This is basically because Hebbian learning builds 
"shortcuts" . If two neurons are not wired (thus far apart from the synaptic graph point of view) but if they are "synchronized" , (e.g. 
i is active at time 4 + 1 whenever j is active at time t) then the learning rule will construct a synapse between them. Note that, as 
discussed in the section IV and in the next section, 2 neurons can be synchronized even i they are not wired, by the mere effect of the 
non linear dynamics (see section VI. F for a discussion of this aspect in chaotic networks). 
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level of activity of active neurons (and inhibit more and more the silent ones), the derivatives of the transfer function 
of the neurons has a tendency (on average) to become smaller. Thus, the entry of A{x{t)) become smaller on average. 
If we (roughly!) replace the condition gop{J'A{x*)) = 1 by gop{J < A(x) >) = 1 and if we neglect the modifications 
of the Jij's induced by learning, one sees that the g value to destabilize the network increases while learning is 
performed. Thus, the effective motion induced by learning in the parameter space corresponds to get closer and closer 
from the destabilization manifold, with an eventual crossing when learning is to long. Finally, since for large N, the 
destabilization manifold and the edge of chaos are very close one concludes that learning lead the system closer and 
closer from the edge of chaos. 

What about presentation or removal of a pattern ? The learning rule (122) depends on the activity of the neurons. 
Since the initial presentation of the pattern leads to changes in the distribution of the neuronal local fields Ui, this 
activity is (possibly slightly) modified by a pattern presentation. Prom the parameter space point of view the pattern 
presentation corresponds to a translation in the subspace of thresholds, in the direction of the vector ^. The whole 
learning phase is conditioned by the presence of the pattern. It has the effect of increasing the numbers of saturated, 
Xi = lov silent, Xi = neurons. The global effect is similar to having an effective threshold whose value grows during 
the learning phase, leading to the observed dynamical collapse. Removing the pattern has in general the effect of 
reducing the width of the distribution of neural local fields and the number of saturated/silent neuron decreases. If 
the system is close to the edge of chaos (this happens when we stop learning slightly after the reduction of chaos to 
a periodic or quasi periodic attractor) this can induce the drastic dynamical change observed. Thus, this scenario 
lead us to conclude that the learning dynamics leads the system "to the border of the chaos", in a state where it is 
sensitive to the learned pattern (i.e. a translation in the direction of the pattern (presentation, resp. removal) induce 
the crossing of the border of the chaos). In some sense, the network has adapted itself to the pattern, via learning. 
This has an interesting echo with the adaptation condition (78) of the section IV. In particular, there should exist 
transversality conditions ensuring that the presentation/removal of the pattern leads to a "transverse crossing of the 
edge of chaos" . 

This discussion gives us interesting hints but is not entirely satisfactory. Firstly, as already said, we don't have 
a real theory to validate this scenario. Also, we didn't discuss the effect of presenting another pattern, after the 
learning phase. More generally all the discussion above dealt with a specific example of a specific rule. What about 
the genericity of this result ? What about its robustness ? What happens if one changes the learning rule ? 

Actually, the rules (122) is rather rough and not really robust. It has been introduced as a straightforward 
implementation of the recipes in section III.C providing an interesting pedagogical example. However, to have robust 
effects one needs to consider more elaborated rules. Systematic investigations have been performed in (134), (51), (50). 
Various learning rules have been proposed, having the general form (58) 

4 = J,, + -r,,. »,j = l..iV; (123) 

where Fjj may either depend on the value of the "instantaneous" pre- {t) and post- synaptic {t + 1) neuron or on 
averages such has F^ = ruirnj or F^ = Cij(l), where Cij(l) is the time 1 correlation function between j and i. 
In the case where Fjj depends on average values, one has to consider two coupled dynamical systems. A fast one 
corresponding to the neurons evolution and a slow one corresponding to the evolution of the J^^'s. In the joint 
evolution one has then to wait that the fast neurons dynamics settle onto its attractor before performing one learning 
step. 

The main observations above remain (50; 132). Moreover, it is possible to improve the learning rule so that 
the response of the system to the pattern in terms of chaos reduction is selective and specific. Presenting another, 
completely distinct stimulus, does not lead to a dynamical reduction. However, a weakly noisy version of the 
stimulus has this effect. Henceforth, this mechanism is robust to noise. It is also possible to learn several stimuli but 
the storage capacity of the learning rule (123) is weak. More elaborated versions can reinforce the storage capacity (51). 

These results are fascinating for they are the demonstration of an effect similar to Freeman's paradigm (even one 
should take care when drawing biological conclusions from this simple model). To the best of our knowledge this 
model is the first example of a formal neural network exhibiting this effect. However, one may can ask what are the 
potential applications of this. Actually, one may complain that to observe this dynamical reduction one needs to 
somehow "assist" learning since one has to stop it before the dynamics irremediably die. Also, learning left us with an 
association pattern/attractor, but how can we use this ? In fact, the more interesting observations are on one hand 
that the spontaneous dynamics is chaotic and on the other hand that learning a given stimulus leads to a repartition 
of active neurons that depend on the stimulus (138). Chaos allows the spontaneous dynamics of the neural network 
to explore a wide range of "possibilities" each them corresponding to a state of the network, while having neurons 
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selectively responding to stimuli/patterns can be used to perform tasks or make decision. In this sense, the Neural 
network (102) can be used as a first layer of a complex neural architecture. This has been for example used in the 
training of an autonomous robot designed to adapt its motion to a random environment (138). 

This subsection leave us with an interesting conclusion. The hebbian learning rule (123) allows us to store some 
information in the chaotic neural network (102) and this information can be somehow retrieved. But this leads 
to several natural questions : How can a chaotic network store and treat information ? Where is the "learned" 
information stored ? Is there a way to see that such a network has learned something without presenting the pattern 
? The collapse effect is clearly a collective effect, but this does not mean that all neurons play the same role in the 
dynamics ? These questions will not be answered in this paper but in the next subsection we present a new analytical 
tool that may, in the long term, be used to tackle such problems. Recent developments have been recently made in 
this direction in (40). 

F. Influence of a time dependent Input: signal propagation and linear response theory. 

Let us first return to the point raised in section III.B. Since synapses are used to transmit neural fluxes (spikes) from 
a neuron to another one, the existence of synapses between a neuron (A) and another one (B) is implicitly attached 
to a notion of "influence" or causal and directed action. However, as we saw, a neural network is a highly dynamical 
object and its behavior is the result of complex interplays between the neurons dynamics and the synaptic network 
structure. Moreover, the neuron B receives usually synapses from many other neurons, each them being "influenced" 
by many other neurons, possibly acting on A, etc... Thus the actual "influence" or action of A on B has to be 
considered dynamically and in a global sense, by considering A and B not as isolated objects, but, instead, as entities 
embedded in a system with a complex interwoven dynamical evolution. In this context it is easy to imagine examples 
where there is a synapse from A to B but no clear cut influence, or, in the opposite, no synapse and nevertheless an 
effective action. 

Consider indeed the figure 36. Neuron 1 excites neurons 3, but in the same time it excites neuron 5, which inhibits 
neuron 3. What is the effective action of 1 on 3 ? This clearly depends not only on the synaptic weight, but also on 
the state of the neurons 1,3,5. More generally, the spikes or signals emitted by a neuron can follow different paths, 
and its effective influence results from the contribution of all these paths. Actually, one can easily figure out by a 
simple glance at figure 36 that feedback loops (that is closed circuits in the synaptic graph) play an important role. 
However, as pointed out several times in this chapter one has to consider topological aspects (such as the feedback 
circuits) and dynamical aspects. 




FIG. 36 Example of network illustrating the eflfect of feedback loops. 

One way of doing this is to compute cross correllogramms. Indeed, the time correlation function CAB{t) between 
the "state" of A and the state of B incorporates the dynamical evolution and the effective effects due to the neural 
network as a whole. However, correlations functions do not really provide causal information. Indeed, a strong 
correlation between A and B at time t does not tell us if A acts on B or if B acts on A (note in particular that 

CAB{t) = CBA{-t)). 

Another way to measure a causal action consists in exciting neuron A, say with a weak signal, and observe the effects 
on B, e.g. by comparing its evolution with and without the signal applied on A. We shall give later on an explicit 
way to do this. Nevertheless, there is a common wisdom, coming from non equilibrium statistical mechanics, stating 
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that the response of _B to a weak perturbation on A (hnear response), if it exists should actuaUy be a correlation 
function. This is the celebrated fluctuation-dissipation theorem (FDT). We shall however see below that the FDT 
may not hold in simple neural networks models, due in particular to saturation effects in the spike rate emission. 

Finally, a natural choice for an excitatory signal is a periodic signal, with a tunable frequency. Thus, the response 
function, drawn versus frequency, provides similar information as the complex susceptibility in physics. In particular, 
peaks in the susceptibility corresponds to resonances, that is a response of maximal amplitude. We shall see below 
how these resonances can be used to provide an effective, frequency dependent notion of network structure. We shall 
also see how they incorporate non linear effects in the dynamics even though they are obtained in the context of 
linear response theory. 

With these ideas in mind consider the model (102) in the chaotic regime and assume that we superimpose upon the 
state Uj (t) of the node j a small external signal (t) . How does this signal propagate inside the network ? Because of 
the sigmoidal shape of the transfer functions the answer depends crucially, not only on the connectivity of the network, 
but also on the value of the UkS. Assume, for the moment and for simplicity, that the time-dependent signal ^j{t) 
has variations substantially faster than the variations of uj. Consider then the cases depicted in Fig. 37. In the first 
case (a) the signal ^j{t) is amplified by /, without distortion if ^j{t) is weak enough. In the second case (Fig. 37b), it 
is damped and distorted by the saturation of the sigmoid. More generally, when considering the propagation of this 
signal from the node j to some node i one has to take into account the level of saturation of the nodes encountered 
in the path, but the analysis is complicated by the fact that the nodes have their own dynamical evolution (Fig. 38). 
A mathematical formulation of this is given e.g. in eq. (128) below. This shows once again that the analysis of this 
signal propagation must take into account the topological structure of the graph as well as the nonlinear dynamics. 




FIG. 37 Nonlinear effects induced by a transfer function with a sigmoidal shape on signal transmission. Fig. 37a. Amplification. 
Fig. 37b. Saturation. 



k2,t=2 




FIG. 38 The propagation of a signal along a path in the network depends not only on the weights of the links but also on 
the level of saturation of the nodes that the signal meets. The level of saturation depends on the current state of the node 
(schematically represented as a red point in the figure). This state evolves with time. 
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In this context wo would like to measure the average "influence" of neuron A on neuron B (namely how a weak 
signal applied on A perturbs on average the state of B), including the effects of the non linear dynamics. There 
is a natural notion of average in chaotic systems such as (102) related the so-called Sinai-Ruelle-Bowen measure p 
(SRB) (143) (see appendix) which is obtained as the (weak) limit of the Lebesgue measure ^ under the dynamical 
evolution^^: 

p= lim G>. (124) 

n — *+oo 

In the following we will assume that all Lyapunov exponents are bounded away from zero. Then for each u G suppp, 
where suppp is the support of p, there exists a splitting E^^ ® E^^^ such that E^\ the unstable space, is locally 
tangent to the attractor (the local imstable manifold) and the stable space, is transverse to the attractor (locally 
tangent to the local stable manifold). Let us emphasize that the stable and unstable spaces depend on u (while the 
Lyapunov exponents are fj, almost surely constant). Let us consider a point u on the attractor and make a small 
perturbation 6^- This perturbation can be decomposed as = <5u + where G -Bu"^ and S'^ G E^\ 6^ is 
locally amplified with an exponential rate (given by the largest positive Lyapunov exponent). On the other hand (5* 
is damped with an exponential speed (given by the smallest negative Lyapunov exponent) 

Assume now that we superimpose a signal of weak amplitude upon some of the "membrane potentials" {uk) in such 
a way that the dynamics is still chaotic (with only a tiny variation of the Lyapunov exponents). (This means that the 
method of signal injection is intended to be non invasive). For simplicity, we suppose that the signal does not depend 
on the state of the system, but we can consider this generalisation without difficulty (linear response still applies in 

this case, but the equations (126,127) do not hold anymore). Denote by ^ the vector {^j}^^. The new dynamical 
system is described by the equation: 

u(t+l) = G[u(f)]+^(t) (125) 

The weak signal ^{t) may be viewed as a small perturbation of the trajectories of the unperturbed system (102). At 
each time this perturbation has a decomposition ^{t) = ^^^\t) + ^^"\t) on the local stable and unstable spaces. The 
stable component ^^*^(f) is exponentially damped. The unstable one ^^"^(f) is amplified by the dynamics and then 
scrambled by the nonlinear terms. Consequently, it is impossible to predict the long term effect of signal ^(t) on the 
global dynamics. 

This is true for individual trajectories. However, the situation is substantially different if one considers the average 
effect of the signal, the average being performed with respect to the SRB measure p of the unperturbed system. 

Indeed, as an application of the general theory (131), it has been established in (37), (38) that the average variation 
{t) of the membrane potential Uj under the influence of the signal is given, to the linear order, by: 

oo 

iMt) - Ui{t)) = Y^Y^^ijiaMt - <7 - 1) (126) 

(7=0 j 

We used the shortened notation < > for the average with respect to p. In this expression Xij{(^) are the matrix 
elements of : 

X{a) = I p{du)DGl (127) 

Thus x{<^) is a matrix representing the average value of the iterate a of the Jacobian. Let us note that the fact 
that x(o") stay bounded for (t ^ oo is not a trivial result because -DG^ diverges exponentially with a. The 
convergence of x{'^) has been rigorously shown by Ruelle under the hypothesis of uniform hyperbolicity. It results 
from the exponential correlation decay (mixing) in the unstable directions and on the exponential contraction In our 
framework, this means that, provided that ^(t) is sufficiently small, and for any smooth observable A, the variation 
< A >t — < A > is proportional to ^(<) up to small non linear corrections. In other words, pt is differentiable with 
respect to the perturbation. The derivative is called the linear response. 

It is interesting to note that the response at time 1 is < DG{u) >, namely this is the average value of the 
Jacobian matrix. Thus, at time 1 we have a complete correspondence between the notion of influence discussed in the 



A crucial property is that a SRB measure has a density along the unstable manifolds, but it is singular in the directions transverse to 
the attractor. This feature is at the origin of the distinction between unstable and stable poles of the susceptibility (see below). 
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section V.C and the linear response. Tliis suggests us to eonstruct circuits of influence as we did in the section V.C. 
Unfortunately, this correspondence does not hold for larger times. This is basically because the quantity < DG'^(u) > 
does not obey the chain rule (contrarily to £>G^(u)). Therefore, if j influences i and if i influences k, this does not 
imply that j influences k. In the case of dynamical system (102) one can decompose XijiT) '■ 



The sum holds on each possible paths "fijir), of length r, connecting the neuron fco = j to the neuron = i, in 
T steps. One remarks that each path is weighted by the product of a topological contribution depending only on 
the weight Jjj and a dynamical contribution. Since, in the kind of systems we consider, functions / are sigmoid, 
the weight of a path jijlr) depends crucially on the state of saturation of the neurons ko,...,kr-i at times 
0, ...,T — 1. Especially, if f'{uk,_-^{l — 1)) > 1 a signal is amplified while it is damped if f'{uk,_i{l — 1)) < 1- 
Thus, though a signal has many possibilities for going from j to i in r time steps, some paths may be "better" 
than some others, in the sense that their contribution to Xiji^) is higher. Therefore eq. (128) underlines a key 
point. The analysis of signal transmission in a coupled network of dynamical neurons with non linear transfer func- 
tions requires to consider both the topology of the interaction graph and the nonlinear dynamical regime of the system. 

One can decompose the response function (127) into two terms. The first one is obtained by locally projecting the 
Jacobian matrix on the unstable directions of the tangent space. This term will be named the "unstable" response 
function. It corresponds to linear response of the system to perturbations locally parallel to the local unstable manifold 
(roughly speaking perturbations "parallel to" the attractor). One can show that the linear response associated with 
this type of perturbation is in fact a correlation function, as found in standard fluctuation-dissipation theorems (131). 
Hence, as usual for correlation functions of a chaotic system, it decays exponentially (because of mixing) and the decay 
rates are associated with the poles of its Fourier transform. More precisely, these exponential decay rates correspond 
to the imaginary part of the complex poles of the unstable part of the susceptibility (128). Thus they will be called 
"unstable" poles. More generally, it can be shown that these poles are also the eigenvalues of the operator governing 
the time-evolution of the probability densities (which we denoted above as G*/i), the so-called Perron- Frobenius 
operator (124). Therefore, these poles, whose signatures are visible in the peaks of the modulus of the correlation 
functions, do not depend on the observable, though some residues may accidentally vanish for a given observable. 

The second term is obtained by locally projecting the Jacobian matrix on the stable directions of the tangent 
space. It corresponds to the response to perturbations locally parallel to the local stable manifold (namely transverse 
to the attractor). Therefore, it is exponentially damped by the dynamical contraction. [Note that, according to the 
specific form of the Jacobian matrix, this contraction is, in our case, mainly due to the saturation of the sigmoid 
transfer function] . The corresponding exponential decay rates are given by the complex poles ( "stable" poles) of the 
stable part of the complex susceptibility. But here the poles depend a priori on the observable. One can easily figures 
this out if one decomposes the stable tangent space of a point in the orthogonal basis of Oseledec modes (directions 
associated to each of the negative Lyapunov exponent). The projection of the i-th canonical basis vector on the fc-th 
Oseledec mode depends on i and k. This dependence persists even if one takes an average along the trajectory, as in 



Hence, both stable and unstable terms are exponentially damped, ensuring the convergence of the series (126), 
but for completely different reasons. Moreover, the stable and unstable part of the linear response have drastically 
different properties. As a matter of fact, the stable part is not a correlation function and it does not obey the 
fluctuation- dissipation theorem. In particular, the unstable poles and stable poles are usually distinct. Moreover, the 
stable poles allow to distinguish the neurons in their capacity to transmit a signal. 

The existence of this linear response theory opens up the way to applications involving chaotic networks used as a 
linear filter. Indeed eq. (126) describes a linear system which transforms an input signal ^{t) of small amplitude into 
an output signal (wj(t) — Ui{t)) according to a standard convolution product. In particular, if the external signal is 
chosen as: 




(128) 



(127). 



(129) 



Note that a linear response theory has also been proposed in (153). However, it requires the invariant measure to have a density. This 
is only true for the conditional measure along unstable manifolds. As a matter of fact, this theory does not contain the stable term. 
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(where is the unit vector in direction j), then the response of the system is also harmonic with : 

iUt) - u,{t)) = ex,,(c^)e-'"^(*-i) (130) 
where the frequency-dependent ampHtude: 

oo 

X.,(c.) = ^X.,('T)e^'^" (131) 

(T = 

is cahed the complex susceptibility. In ref.(37) we have conceived and implemented a method to compute Xiji^) 
numerically. The knowledge of the susceptibility matrix is very useful as it enables one to detect resonances, i.e. 
frequencies for which the amplitude response of the system to a periodic input signal is maximum. In fact the 
existence of a linear response implies that Xij{w) is bounded for all w G [0, 27r]. Moreover, in view of eq. (131), it is 
analytic in the complex upper plane. On the other hand, Xij(uj) can have poles within a strip in the lower half plane, 
e.g. in Wo — iA, A > 0. In this case, and if A is small, the amplitude ('i')! exhibits a peak of width A and height 
|Xij(wo)| which can be interpreted in the present context as follows: when unit j (whose state varies chaotically due 
to the global dynamics) is subjected to a small periodic excitation at frequency uq and amplitude e then the average 
response of unit i behaves periodically with same frequency and amplitude e|xij(ti;o)| which is maximal in a frequency 
interval centred about ujq. 

Let us numerically computes the susceptibility x{uj) for real values of ui (see (37), (38) for details) in the following 
example. This is a sparse network where each unit receives connection from exactly K = 4 other units (sparse neural 
networks of type (102) exhibits also chaos via quasi periodicity (57)). The number of units was fixed to A'' = 9. 
The Jij's have been drawn at random according to a Gaussian distribution with mean zero and a variance The 
corresponding network is drawn in Fig. 39. (Note that the corresponding graph is not decomposable). Blue stars 
correspond to inhibitory links and red crosses to excitatory links. In this example the unit 7 is a "hub" in the sense 
that it sends links to almost every units, while 0, 2, 3 or 5 send at most two links. 




FIG. 39 Connectivity matrix (Fig. 39a, on the left) and the corresponding network for the investigated system (Fig. 39b, on 
the right). In Fig. 39b each node is represented by a circle. A filled circle means that there is a link from the corresponding 
node to itself (red: self-excitation, blue: self-inhibition). Inhibitory links are terminated by a vertical bar while excitatory links 
are terminated by an arrow. 

A small constant 0i has been added to each Ui to break down the symmetry u — > — u (i.e. Ui(t) = JijXj{t) + 9i). 
As expected the corresponding dynamics exhibits a transition to chaos by quasi-periodicity. For 5 = 3 the dynamics 
has one positive Lyapunov exponent (Ai = 0.153) and 8 negative Lyapunov exponents (with A2 = —0.427). The 
Lyapunov exponents have been computed with the Eckmann-Ruelle algorithm (58). The chaotic regime is stable to 
small perturbations, as we checked. 

Computing the susceptibility one obtains the curves shown in Fig. 40. Several remarks can be made. First, some 
resonance peaks are rather high {'^ 20) corresponding to an efficient amplification of a signal with suitable frequency. 
It is also clear that the intensity of the resonance has no direct connection with the intensity or the sign of the 
coupling and is mainly due to nonlinear effects. For example, there is no direct connection from to 3 or 5 but 
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FIG. 40 . Modulus of some susceptibilities. Fig. 40a (top). 7 (highly connected unit) excites the units: 1 (excitatory link 
with intensity J17 — 0.007); 2 (no direct link); 3 (excitatory link with intensity J37 = 0.722); 6 (inhibitory link with intensity 
Jq7 = —0.041). Fig. 40b (middle). (weakly connected unit) excites the units: 1 (inhibitory link with intensity Jio = —1.131); 
3,5,8 (no direct link); Fig. 40c (bottom). 5 receives the excitation from the units: (no direct link); 1 (excitatory link with 
intensity J51 — 1.015); 5 (no direct link); 6 (inhibitory link with intensity J56 — —1.312). 



nevertheless these units react strongly to a suitable signal injected at unit 0. 

Let us now compare the Fourier transform of the correlations function Cij{t) for the same pairs (Fig. 41). One 
remarks that these functions exhibit less resonance peaks. This is expected since the Fourier transform of the 
correlation function Cij{t) only contains unstable resonances while the susceptibility contains stable and unstable 
resonances. Note also that the resolution in resonance peaks is quite better in the susceptibility. 

The previous analysis leads then us to propose a notion of "effective", frequency dependent, connectivity based on 
susceptibility curves. For a given frequency w, we plot the modulus of the susceptibility |xij(w)| with a representation 
assigning to each pair i, j a circle whose size is proportional to the modulus. Some examples are represented in Fig. 
42. We clearly see in this figure that changing the frequency changes the effective network. 

For example, with a frequency ui = 0.125 (Fig. 42a), the node 1 has a strong ability to transmit signals towards the 
node 5 (namely the response of this unit is high). On the contrary, nodes 5, 6 and 7 have weak performances in signal 
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FIG. 41 . Modulus ol the correlation functions corresponding to the susceptibilities represented in the Fig. 40a,b,c. 



transmission at this frequency. Moreover, one sees that 7 is a bad sender and a bad receiver. With a frequency 0.57 
the effective network has a rather symmetric structure and basically all units respond to this excitation (however with 
a different amplitude). Also, some units present a strong affinity with some others, at a specific frequency. Obviously, 
one also checks that for frequencies that do not correspond to resonances (such as = 2.33 in Fig. 42f) the response 
is essentially inexistent whatever the pair. Finally, this figure shows that it is possible to excite any unit from any 
other one in such a way that this unit (and possibly a few other but not all the other units) have a maximal response. 

All these effects are due to a combination of topology and dynamics and they cannot be read in the connectivity 
matrix J^. 



G. Conclusion 

This section was devoted to the analysis of some recurrent neural networks, which are a particularly prominent 
example in this field. We have analyzed in some details the collective dynamics and exhibit several important effects 
revealing the richness and complexity of the emergent dynamics. Indeed, as noted in the begin of the section, the 
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FIG. 42 Effective connectivity for : Fig 42a (top left) uj = 0.125; Fig 42b (top right) uj = 0.57; Fig 42c (middle left) uj = O.i 
Fig 42d (middle right) lo = 1.0; Fig 42e (bottom left) co = 2.3; Fig 42e (bottom right) lo = 3.14. 



dynamics of the uncoupled neurons is rather poor. This justifies somewhat the claim, made in the introduction, 
that one can make rather drastic simplifications in the description of the neurons of a coupled system, and still get 
a complex and relevant model. However, one must be cautious. Removing some characteristics and still get an 
interesting behavior does not mean that the removed characteristics are irrelevant "details" . Actually, the models 
presented here are quite simplistic as "brain" models. To our opinion, there main interest is to provide "benchmarks" 
for developing and testing tools that one may use, later on, to analyze more realistic models. 



VII. CONCLUSION 

In this chapter, we have provided examples suggesting that the mathematical analysis of neural networks dynamics 
can be pushed relatively far, in some simple models. However, a remaining question is: can we perform the same kind 
of analysis for neural networks closer to biological systems ? At the actual stage of research the techniques of cerebral 
imagery; brain analysis and neurophysiology allows to go relatively deep in the structure and dynamics of cerebral 
areas, but it allows also to make an explicit cartography of the nervous system of primary animals such as worms 
(e.g. Caenorhabditis elegans (31)). Thus, it is in principle possible to write the explicit dynamical system accounting 
for the evolution of small area containing a relatively small number of neurons (^ 100 — 1000). However, the detailed 
analysis of these equations is still intractable. This is also, in some cases, useless. Indeed, often these area exhibit a 
relatively simple collective behavior. It is thus profitable to define a phenomenological model, described by a small 
set of differential equations and a few parameters that one can adjust to fit experimental results. 

A prominent example concerns the cortical columns implicated in vision. A cortical column is a population of 
pyramidal cells receiving excitatory and inhibitory inputs from others cells in the same column but also excitatory 
inputs coming from other columns, close or distant. A celebrated dynamical model has been proposed by Lopes 
Da Silva (108) and Jansen & Rit (98) describing the activity of cortical columns. A mathematical analysis of the 
bifurcations exhibited by this model has been performed in (77). It shows the existence of oscillations generated by a 
Hopf bifurcation, induced by the variation of a parameter modeling the frequency of stimuli emitted by an external 
source. The value of the oscillations frequency is about 10 Hz corresponding to the a rhythm. One can also exhibit 
spikes emission, looking very much like epileptic activity, and related to a saddle-node bifurcations on a cycle. In 
the vicinity of this bifurcation, an excitation with an external stimulus with a specific frequency induces a spike 
train emission. The techniques used by the Authors combine standard results from dynamical systems theory and 
numerical analysis, in the spirit of the analysis presented in the section II. C. Note that rhythm (5, 0,/3, 7 have also 
been numerically exhibited in (49), when varying the excitatory-inhibitory effects in biologically realistic ranges. 

This example shows that it is indeed possible to analyse neural networks closer to biology, possibly after some 
simplification of the initial system. More details about vision and cortical columns are given in chapter IV. 
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We have made a trip in the world of Neural Network dynamics, following the path represented in Table 1. As 
said in the beginning many examples, models, etc . . . have been omitted. However, we have tried to give an outlook 
of the various methods available for the study of the dynamics. This excursion has also shown that, when going 
from a level of complexity (one neuron dynamics) to another level (collective dynamics), it might be fruitful to adopt 
different perspectives (accurate description of a neuron versus emergent behavior of "simplified" neurons) and different 
(but complementary) methods (dynamical system theory versus probability theory and statistical physics). It also 
shows us the necessity to develop accurate tools to handle neuronal dynamics (this is well known and not new) and 
the possibility to do this by combining existing theories and numerical analysis. This is a formidable task but the 
byproducts are on one hand a better understanding of neuronal dynamics and on the other hand a possible insight in 
other fields. 
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VIII. APPENDIX 

This appendix is mainly devoted to non-specialists. It gives a brief summary of the concepts and techniques in 
dynamical systems theory used in this chapter. Our main references are (12), (13), (78), (100), (129). 

A. Elementary notions in dynamical systems theory. 

1. Basic definitions. 

The dynamical systems studied in this chapter are either defined by a (finite) set of differential equations : 

^=H(X;A) (132) 

or a set of recurrences^*: 

X(t + 1) = F(X; A) (133) 

where X e A1, Al being a compact set in R''^, where A'' is the number of degree of freedom and X denotes the vector 
{a;^}^^. The vector fields H (resp. the recurrence F) in eq. (132), (resp. (133)) do not depend explicitly on time. 
The corresponding dynamical system is then called autonomous. We mainly deal with the autonomous case in the 
paper and in this appendix. \ E £\ d refers to a set of p (real) parameters on which the system depends. This 
might be an external current applied to a neuron, an external input submitted to an assembly of neurons, the set of 

synaptic weights, etc Therefore, A can have a large (though finite) dimension. It can also be deterministic or 

random. The last case requires however combinations of techniques from dynamical system theory and probability 
theory. An example is developed in section VI. 

We assume that H, F are smooth (at least C^) functions of X, A. In the continuous time case (132) the Cauchy 
theorem ensures the local unicity of the solutions provided that H is a Lipschitz function. Namely, if X G , there is 
a time interval ] — c, c[ and a neighborhood W 9 X such that there is a unique solution of (132), X(i) GU,t g] — c, c[ 
and such that X(0) = X. Moreover, when A4 is compact, the solutions extend to t G [— oo,-|-oo[ (43). Denote by 

X =^ {'K{t)}f^ the (forward) orbit or trajectory such that X(G) = X and by x~ {'X.{t)}'^^_^ the backward 
trajectory. The unicity of trajectories implies that two trajectories cannot cross (though they can accumulate on the 
same set, as shown below). Also, the equations (132) have the meaning that any trajectory is locally tangent to the 
vector field H. In low dimensional cases (namely N < 3) this is helpful to draw a qualitative sketch of the main 
dynamical system features (phase portrait), without any computation (sec for example the sections II.B.2,II.D). 

In the case of the recurrence (133) the forward trajectory is simply constructed by iterating the map F. Therefore 
it is always defined (provided that the initial condition is in the domain of definition of F). The backward trajectory 
is uniquely defined only if F in invertible. In the sequel we shall assume that F is a diffeomorphism. For the 
dynamical system (132) one can prove the existence of a one parameter family of diffcomorphisms (or flow), such 
that (/>° = id, ^* o 0** = (^*+* and X(t) = (/>*(X). In the sequel, we shall use the notation X(t) = f*(X) for both dy- 
namical systems (132), (133). Consequently, f will refer to the flow in the case (132) and to the map F in the case (133). 

The dynamical systems (132), (133) may exhibit a wide variety of dynamics, from very simple (rest state attracting 
all trajectories), to complex (chaotic behavior) and even more complex (coexistence of many chaotic attractors, etc 
...). Consequently, in most cases the explicit solution of (132), (133) are not known. The current philosophy in 
dynamical systems theory, initiated by H. Poincare (122), is that finding a general solution is not only impossible, 
but also useless. Indeed, in many cases, a qualitative study of the dynamical system is enough to extract quite a 
large amount of informations which often allows us to capture the main features of the dynamics. In particular, one 
can extract characteristic ensembles such has attractors, repellors, periodic orbits, etc . . . , which contains the main 
informations one needs. In many cases, one is indeed interested in the asymptotic behavior of the forward orbits. The 
ui-limit set of X is the set of accumulation points of the forward trajectory X(f). The w-limit set of f is the union of 
the w-limit sets for all X e A^. It contains in particular the attractors of the dynamics (see the definition below). The 
same notion (a-limit) set can be defined for the backward trajectory when it is defined. A more general and related 
notion is the non wandering set. This is the set of points X such that for any open neighborhood U B H- there is a 



This implies that we do not consider the case of Neural Networks with sequential dynamics. 
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time to > such that {''"{U) HU ^. This set contains the main elements of the dynamical system such as the lo 
limit set. 

The w-limit set lay have a quite complex structure. However, it contains in general some characteristic objects such 
as fixed points, or periodic orbits. X* is a fixed point if its orbit consists of X* only. In other words, H(X) = [resp. 
F(x) = X]. a is a periodic point if there is some t > such that f*a = a. The lower bomid of such t is the period of 
a, T{a). The set T = {f*a ;0 <t < T{a)} is called a periodic orbit or a closed orbit. For a discrete time dynamical 
system it is a finite set; for continuous time it is continuously infinite. 

2. Fixed points and linear analysis. 

The first step of the analysis of (132), (133) is to seek for equilibria or fixed point. A fixed point is stable iff for 
any neighborhood U 9 X*, there exists a neighborhood Ui C U such that VXq E Ui, Vt > 0, f*(Xo) G U. X* 
is asymptotically stable if there exists an open neighborhood Ui such that VXq G Ui, f*(Xo) X* as t ^ oo. 
Asymptotically stable fixed points are called sink. A stable fixed point which is not asymptotically stable is called a 
center (see Fig. 43). A well known example is the stable equilibrium position of the imdamped pendulum. A fixed 
point is unstable if it is not stable. Note that the notion of stability is a local notion. Among the various kind of 




sink center saddle source 



FIG. 43 Various kind of fixed points. 

fixed points, the stability of hyperbolic fixed points can be analyzed by linearization about X*. Indeed, call DHx* 
(resp. DFx*) the Jacobian matrix of H (resp. F) at X*. Since the coefficients of this matrix are real, the eigenvalues 
are either real or complex conjugate. Call S*!? [^] the spectrum of a matrix A. One decomposes Sp[DHx*] [resp. 
Sp[DFx'] ]into three parts : the stable eigenvalues are such that 3?(A) < [resp. |A| < 1]; the neutral eigenvalues 
are such that SR(A) = [resp. |A| = 1] and the unstable eigenvalues are such that §?(A) > [resp. |A| > 1]. Moreover, 
the Jacobian matrix can be reduced to a diagonal (or more generally to a Jordan normal form) in a basis vi , . . . v^v 
corresponding to the (generalized) eigenvectors. The stable space f (X*) is the subspace of IR^ generated by the 
eigenvectors corresponding to the stable eigenvalues. In the same way one defines the central space £^{X.*) and the 
unstable space 5"(X*). 

Then X* is an hyperbolic fixed point of (132) if there is no neutral eigenvalues (resp. 5'^(X*) = 0). X* is linearly 
stable if additionally f "(X*) = (namely all eigenvalues are stable). A linearly stable equilibrium is asymptotically 
stable and the rate of convergence is given by the largest real part of the eigenvalues in the case (132) (continuous 
time), and by the largest modulus of the eigenvalues in the case (133) (discrete time). Unstable hyperbolic fixed points 
are divided into saddle points (there are stable and unstable eigenvalues) and sources (all eigenvalues are unstable) 
(see Fig. 43). Hyperbolic fixed point have the following important properties'^. 

1. Hartman-Grobman linearization theorem. If X* is hyperbolic then there exists an homeomorphism h preserving 
the sense of orbits, locally mapping the orbits of the flow of (132) (resp. the map (133)) to the orbits of the 
linear flow eJ^^^' (resp. the linear map DF^,). The Hartman-Grobman theorem implies that the dynamics 
near an hyperbolic fixed is essentially equivalent (up to a smooth variable change) to a linear system (for a nice 
application to Neural Networks see section IV.B). 

2. Invariant manifolds. Let U he a, neighborhood of X*. If X* is hyperbolic then there exists local stable and 
unstable manifolds: 



Note that the notion of hyperbohcity extends to moving points (see section VIII. C. 2) and that the result below can be generalized (see 
(100)) 



70 



Wil^iX*) ^ {yeU |f*(y) -^X* ast^oo and f*(y) eU,\ft> O} (134) 
W,'^,^(X*) ^ {yeU |f*(y) -> X* as i -> -cx) and f*(y) eU, yt< O} (135) 

respectively with the same dimension ns , n„ as the eigenspace E^, , E^, of the Hnearized system, respectively 
locally tangent to Ej^,,E^, at X*, as smooth as the function H (resp. F) and dynamically invariant. Moreover 
the angle between E^, , E^, is bounded away from zero. The local stable and unstable manifold have global 
analogues : 

W'*(X*) = Ut<oWf„,(X*) = {yeM |f*(y) ^ X* as t ^ oo} (136) 
W"(X*) - Ui<o>Vr„,(X*) = {y e M |f*(y) ^ X* as t ^ -ex.} (137) 




FIG. 44 Fig. 44 a. Hartman-Grobman theorem. Fig. 44 b. Local stable and unstable manifolds. 

Stable and unstable manifolds may intersect in homoclinic (Fig. 45a, b) or heteroclinic intersections (Fig. 45c) . This 
has important consequences. In particular, the global stable and unstable manifolds of a fixed point may have strong 
influence to the global dynamics, as in the case of transverse homoclinic intersections for maps ((14; 78)). 




FIG. 45 Fig. 45 a. Homoclinic intersection (continuous time system). 45 b. Transverse homoclinic intersection (discrete time 
systems). Fig. 45 c. Heteroclinic intersection. 

The notion of fixed point is related to a more general notion called "convergence". Following Hirsch (87) we say 
that a dynamical system is 

• Convergent: if all trajectories converge to equilibria. 

• Globally convergent or asymptotically stable: if all trajectories converge to a unique equilibrium. 
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3. Lyapunov functions. 

A Lyapunov function V is a differentiable^° function which decreases along the trajectories and is bounded from 
below. In dissipative mechanical systems, the energy is a Lyapunov function. This notion is useful to locate fixed 
points (they arc cxtrcnia of V) and to analyze their stability. Indeed if ^ < (resp. ^ < 0) in the neighborhood of 
some fixed point X* then X* is stable (resp. asymptotically stable.) More generally the Lasalle invariance principle 
(106) asserts that the w-limit set of any point X is included in the largest invariant set where V is a constant. An 
important corollary is that if V is a strict Lyapunov function (^ < 0) on a compact set M. then the equilibria are 
isolated, and the system is convergent. Lyapunov functions are used in section V.B. 



B. Bifurcations. 

The dynamical systems (132), (133) depend smoothly on a set of parameters A G £\. When varying these parameters 
one modifies the dynamics. On open domains of parameters the changes are essentially quantitative, namely a variables 
change maps the initial system to the modified one. One says that two flows (or maps) f ,f' are topologically equivalent 
if there is an homeomorphism mapping the orbits of f to the orbits of f and preserving the ordering of points along 
the orbits. Two topologically equivalent dynamical systems have therefore the same phase portrait (but quantitative 
characteristics such as the convergence rate to a fixed point may differ). A dynamical system f is structurally stable 
if any sufficiently close f ■^^ is topologically conjugated to f . 

There exists in general a (closed) set of parameter values where the corresponding dynamical system is not 
structurally stable. At these points, called bifurcations points, the dynamics changes qualitatively. The codimension 
of the bifurcation is the number of independent parameters one has to adjust in order to obtain the bifurcation. 
In this section we focus on bifurcations oc;curring on fixed points. Moreover, we only consider the case where at 
most two independent parameters are varying. This is indeed the only cases where a complete classification of fixed 
bifurcations is known (78). 

Assume therefore that X* is a fixed point, namely this is the zero of some function G(X; A) (G(X; A) = -ff (X; A) 
in the continuous time case, and G(X; A) = F(X; A) — X in the discrete time case). When varying A the implicit 
function theorem guarantees that X* moves along a regular curve X*(A) provided that DG(X; A) is invertible. This 
also implies that the eigenvalues of -DG(X;A) are moving continuously. Note that, since Z)G(X;A) is real, the 
eigenvalues are either real or complex conjugated. Then, at some parameter values, some eigenvalues can intersect 
the real axis (resp. the unit circle in the discrete time case). There are two possibilities. Either they cross at the 
origin (resp. at 1). In this case the implicit function no more applies and several branches of solutions of the equation 
G(X; A) = appear or disappear (see Fig. 46,48). Or they cross at imaginary values. This induces in general a 
change of stability for X*(A) and the appearance or disappearance of a limit cycle (see Fig. 49). 

The initial dynamical system has N degree of freedom. However, at the the bifurcation point, say Ac, one expects 
that the only relevant information is contained in the eigendirections corresponding to the crossing eigenvalues. This 
leads to a general method called the central manifold reduction. Let S^(X*) be the central space (it is non zero at 
the bifurcation point), Uc = dim{E''{'X*)) the number of crossing eigenvalues and call £'''(X*) = i?*(X*) _E"(X*). 
Then the central manifold theorem (30) states that there is a function ?i(X;A) : E'^ x S\ ^ E'^{'K*) such that 
W(X*, Ao) = 0, D^n{X*, Ao) = and such that the manifold: 

W%X) = {X + n{X; A) I X e } 

contains X* and is tangent to E^, at this point. Moreover VV'^(A) is locally invariant for A sufficiently small and 
bounded. This means that there is an open neighborhood U of X* such that if X(0) G W{X)nU then X(i) G yV"=(A) 
as long as X(t) G U. Finally >V'^(A) is locally attractive if i?"(X*) = 0. Therefore, in this case, all solutions staying in 
U tend exponentially fast to some trajectory on H'^(A). W^(A) is called the center manifold (though it is not unique). 



Note that the continuity is sufficient in the definition. However, difFcrentiability allows us to replace the condition ^ < by ( VV, H) < 0, 
where <, > is a scalar product in R^. This means that the trajectories cross the level curves of V "inward". Note that, reciprocally, if 
there exists a metric such that < VV,H >< 0, V is a Lyapunov function for the corresponding dynamical system. This allows one to 
show the convergence of some dynamical systems under quite general conditions (see section V.B and (24).) 
See (129) for a definition of a topology in a space of flows. 
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It is then possible to locally reduce the dynamics (60) to the dynamics on W"^(A) by projection. Denote by 
n= : IR^ E" the projection onto E", by n'' : IR^ ~> E''(X*) the projection onto E''{X*), and set X=(t) = n=X(i). 
Then if X{t) is a solution of (132) such that X{t) € W%X) nU,t>0 one has X%t) = X{t) +niX{t);X), namely, 
by a suitable (local) variable change one can write down a smaller dynamically system leaving on W^X) and 
characterizing the relevant part of the dynamics about X* . 

It is then possible to further reduce the dynamics by removing some non linear terms with the appropriate variable 
changes. Actually, one cannot remove all the non linear terms in this way (otherwise the dynamical system is basically 
a linear system). Only the non linear terms satisfying non resonant conditions (see (13; 78) for details) can be removed. 
Finally, one ends with a set of canonical equations called a normal form. In some sense, the normal form reduction 
for a dynamical system is a generalization of the diagonalisation for a matrix. There are uncountably infinitely many 
matrices in IR^ but many matrices have the same diagonal (or Jordan) form. This means that they are equivalent, up 
to a basis change, and the canonical form of their equivalence class is the Jordan form. In the same way, an infinite 
number of dynamical system undergoing a bifurcation at a fixed point can be represented under a canonical form or 
normal form. 

It is remarkable that the different possible codimension one and two bifurcations are in fact only a few. Moreover, 
it is possible to write down general conditions on the dynamical system, called transversality conditions, allowing to 
characterize the type of bifurcation occurring. We now briefly describe these bifurcations. 



1. Codimension one bifurcations. 

In this section, we assume that X* = is a fixed point, and that A is one dimensional parameter. We review now 
the bifurcations arising generically in this case. We denote by Ao the parameter value where the bifurcation arises. 
We first consider the continuous case, and then the discrete time one. 

• Saddle-node bifurcation. The transversality conditions, when written in a great generality, are rather ab- 
stract. However, it is easy to understand them by taking a one dimensional example. Consider indeed the 
system x — f{x; A) such that a; = is a fixed point, and Aq = is a bifurcation point. Performing a Taylor 
expansion about (0; Ao) gives: 

fix; A) = /oo + fiox + /oiA + fnxX + /aox^ + . . . (138) 

Since we want to characterize the dynamical system in a neighborhood of (0, 0) it is natural to consider the 
lowest order terms. Since is fixed point /oo = 0. Moreover, Ao = is a bifurcation point where fj(0; 0) = 
which implies /o,i = 0. If we ask now that the linear term in A does not vanish we get the first transversality 
condition for the saddle- node bifurcation (in one dimension): /qi = f{(0;0) ^ 0. At the bifurcation point the 
implicit function theorem does not apply and two branches of equilibria emerge (or disappear), with a vertical 

tangent (see Fig. 46). The second transversality condition, /20 = §^(0; 0) ^ ensures that these curves have a 
quadratic tangency at (0;0). One can then show that the term /qi and all higher order terms in the expansion 
(138) can be removed. One can also use variable changes which allows us to eliminate /lo, /20- We finally end 
up with the normal form for the saddle-bifurcation in continuous time case. 

x = X-x'^ (139) 
The corresponding bifurcation diagram is drawn in Fig. 46. 

In a A'' dimensional dynamical system the previous discussion generalizes as follows. Assume that : 

(SNl) £>H(X*; Ao) has a simple eigenvalue with right eigenvector v and left eigenvector w. It has also k stable 
eigenvalues and N — k — 1 unstable eigenvalues (counting multiplicity). 

(SN2) 

w.^(X*,Ao)^0 



(SN3) 

' dxjdxk 



cl (X*, Ao)t^jt^fc 7^ 



where we used the Einstein convention (sum over repeated indexes). 
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FIG. 46 Saddle node bifurcation. 



then the normal form of the bifurcation is (139). Namely the dynamical system behaves like eq. (139) in the 
direction of the zero eigenvector, with hyperbolic behavior in the complementary directions. This bifurcation is 
in some sense the "most" generic since the set of dynamical systems which satisfy the transversality conditions 
(SN1),(SN2) is open and dense in the space of C°° one parameter families of vector fields with an equilibrium 
with a zero eigenvalue. 

For discrete time dynamical systems, the normal form writes : 

x{t + l) ^ x{t) + (140) 



• Transcritical bifurcation. Assume now that X* is a fixed point before and after the bifurcation. This implies 
that /oi must be zero and the corresponding transversality (TSN2) condition cannot hold. If we replace it by 
the condition /n we obtain a normal form. 

x^ Xx-x'^ (141) 

whose bifurcation diagram is depicted in Fig. 47. Two fixed points coexist and they exchange their stability at 
the bifurcation point. 

A 

i 

f 




FIG. 47 Transcritical bifurcation. 



In the general case one has to replace the transversality condition TSN2 by: 



(TT2) 



For discrete time dynamical systems, the normal form writes : 

x{t + 1) = x{t) + Xx-x'^ (142) 



• Pitchfork bifurcation. In some cases one has particular symmetries in the dynamical system. A particularly 
prominent example corresponds to the symmetry X — s- —X. Returning to our one dimensional example we see 
that and /2,o has to be zero. We must then consider higher order terms. It is clear that the first remaining 
non linear term is fs^sx^ (and the next one is /s.s). All other terms of order < 3 vanish. Note that according to 
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the sign of /s^s one has supercritical (/a, 3 > 0) or subcritical bifurcation (/a, 3 < 0). In this last case one has to 
take into account the term /s^s in order to "saturate the instabiUty" . The normal form for the (supercritical) 
pitchfork bifurcation is 

x^\x-x^ (143) 
The corresponding bifurcation diagram is drawn in Fig. 48. 




X 



FIG. 48 Pitchfork bifurcation. Fig. 48 a. Supercritical. Fig. 48 b. Subcritical 



• Hopf bifurcation. Assume now that there is a pair of complex conjugate eigenvalues of the Jacobian matrix 
crossing the imaginary axis [resp. the unit circle] at the bifurcation point. Note that this requires that the 
dynamical system has at least a dimension 2. Having eigenvalues with an imaginary part implies that the 
trajectories are locally oscillating around the fixed point. When the eigenvalues cross from the left to the right 
the oscillations are exponentially damped before the bifurcation point, and they are exponentially amplified after 
the bifurcation (see Fig. 49). The exponential amplification is obviously local. When moving away from the 
fixed point the nonlinearities saturate the instability the trajectories converge to a limit cycle. This corresponds 
to a Hopf bifurcation. 

More generally, (78; 112) suppose that the dynamical system (132) as an equilibrium (X*; Aq) such that 

(THl) Z)H(X*; Ao) has a simple pair of pure imaginary eigenvalues ±iuj and no other eigenvalues with zero real 
parts. 

Then there is a smooth curve of equilibria (X(A); A) with X(Ao) = X* . The eigenvalues /i(A), fl(\) vary smoothly 
with A. If, moreover 

(TH2) 

then the normal form is: 

z = jz + az'^z + 0{\z\^) (144) 

where 2 is a complex variable corresponding to the reduction on the central manifold and 7 = (A — Aq) + iuj. 
Note that, written in polar coordinates (r, 6) the equations (144) are: 



f = [d{\ - Xa) + ar^] r (145) 
uj = oj + c{X ~ Xo) + br^ (146) 

where a, 6, c are coefficients depending on the vector field (note that a can be positive or negative, corresponding 
to supercritical or subcritical bifurcation). We remark that the amplitude of the limit cycle increases like the 
square root of the difference A — Aq and that the frequency depends on the parameter and on the amplitude. 
Note also that, near the bifurcation point, the frequency is non zero (see section II. B. 2). The discrete time case 
is substancially more complicated, with specific cases corresponding to strong resonances. A detailed analysis 
can be found in (13). For a discussion in the context of neural networks see also (35) 
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FIG. 49 Hopf bifurcation. 



2. Codimension two bifurcations. 



We describe now one local codimension 2 bifurcations, the Bogdanov-Takens bifurcation. (We only focus on the 
examples found in this chapter). A complete description can be found in (78). Basically, codimension two bifurcations 
may arise either if additional degeneracies in the non linear terms of the previous bifurcations arise, or if the linear 
part of the vector field (the map) is doubly degenerate. In this last case, the linear part for flows takes the form 



1 




• Bogdanov-Takens bifurcation. 

normal form is: 




/ 



V 



-LOl 





CJ2 






-UJ2 





(147) 



The Bogdanov-Takens corresponds to the first situation in eq. (147). The 



X = y 

y = \i + \2y 



(148) 



axy 



where cr = ±1. In the sequel we shall consider the case a ~ 1. The second case can easily be obtained from 
the first one by the substitution t — > —t] y — > —y. It is easy to show that Hopf bifurcation occurs on the curve 
A2 = V— Ai (hence for Ai < 0) while saddle-node bifurcations occur on Ai = 0; A2 ^ 0. The complete bifurcation 
diagram is represented in Fig. 50. 



C. Chaotic motion. 

1. Attractor 

In the previous sections we have seen several example of topological objects where the dynamics converges asymp- 
totically: asymptotically stable fixed points and stable limit cycle are examples of attractor. But attractor can be 
quite a bit more complicated objects. Though there are many different (and non equivalent) definitions of attractors 
(45; 58; 78; 115; 154), they basically combine a notion of attractivity and indecompos ability. Here is one definition 
(100). An attractor is a compact set A such that: 

• (i) Attractivity. There exists an open set U A and a time to such that F*°(W) C U and A = n^oF*(W) 

• (ii) Indecomposability. VX, y ^ A and Ve > there is a chain X = Xq, Xi, . . . X„ = y and a sequence of times 
ti,t2, . . .tri > 1 such that the distance between F*'(Xi_i) and Xi is < e. 

Note that (i) implies the dynamical invariance of A. 

Attractors can have a simple topological structure (fixed points, cycles, tori) or a complex one (strange attrac- 
tors). Though there are several definitions of the "strangeness" of an attractor, there is a general consensus about 
the necessity to have initial condition sensitivity. This notion is in fact related to a more general notions called 
hyperbolicity. 
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FIG. 50 Bogdanov-Takens bifurcation diagram. 



2. Hyperbolic dynamical systems. 

A dynamical system is uniformly hyperbolic if there exists < A < 1 < /i and a constant C such that: 

• (i) There exists two subspaces £^(X),£"(X) respectively called stable and unstable^ forming an invariant de- 
composition of the tangent space at X: Tx = f "(X) f "(X) et Di^E'^iX) = £''(f*(X)) (rcsp. Di^E''{X) = 
f "(f*(X))), Vt > 0, and such that the angle between the two subspaces is bounded away from 0. 

• (ii) Z^fx is contracting on £''(X): If v is a vector in £^{K.) : 

pfxvl! < CA*||v||, > 

• (iii) Dfx is expanding on f "(X): If v is a vector in £"(X) : 

IlL'fx'vll < C/i-*||v||, Vt > 

(Note that the constant C in the definition is independent of X. More generally (non uniform case) this constant 
depends on X. 

(Uniformly) hyperbolic dynamical systems have several remarkable properties (see (100) for a wide description): 
existence of smooth local stable and unstable manifolds locally tangent to the spaces £^(X) (resp. f "(X)); shadowing 
lemma; density of periodic unstable orbits leading to trace formulas; local product structure allowing the construction 
of Markov partition used in symboling coding; structural stability; etc .... 

The existence of an unstable direction implies initial conditions sensitivity while the existence of contracting di- 
rections corresponds to asymptotic convergence onto an attractor. Basically, a strange attractor is composed by the 
closure of the union of the unstable manifolds. A perturbation "parallel" to the attractor (locally tangent to the un- 
stable space) is locally expanded at exponential speed (initial condition sensitivity) while a perturbation transverse to 
the attractor (locally tangent to the stable space) is asymptotically damped. Parallel and transverse time dependent 
perturbations induce drastically different effects on the dynamics (see section VI. F) having interesting interpretation 
in the context of Neural Networks. 

3. Statistical approach and ergodic theory. 

In chaotic systems, it is often useful to replace the study of individual trajectories by a statistical analysis of the 
evolution in the phase space. The natural object is then a probability measure. Actually, there is a close relationship 
between the (physically relevant) probability measures and the notion of state in statistical physics. The initial 

probability /i^*^^ corresponds to randomly selecting the initial conditions and the probability at time t, 
is the result of the action of the dynamics F on It is given by: 

[B] = ^(0) [{X I F*(X) e B]] = [F-*(S)] 
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(149) 



(The corresponding notion in statistical physics is the notion of phase). 

Among all invariant measures the ergodic measures play an important role (they correspond to pure phases in 
statistical physics). There are several equivalent definitions but the most known is certainly the identity between time 
average and ensemble average. A measure /u is ergodic if for /x almost every initial condition X : 



where </> is a function in L^{d^,). 

The definition (150) is unfortunately rather poor since one can show that a dynamical system in a compact space 
has often infinitely many such measures (58; 100). A more useful notion is the Sinai-Ruelle-Bowen (SRB) measure. A 
measure /x is a SRB measure (or natural, or physical measure) if the property (150) holds for a set of positive Lehesgue 
measure (131; 152) of initial conditions. This means basically that the time average and the ensemble average are 
equal for typical initial conditions. Sinai, Ruelle and Bowen have shown that the SRB measure is a "Gibbs like" 
measure: it has an exponential form, although the term in the exponential is not the Hamiltonian encountered in 
statistical mechanics but a dynamically relevant quantity. Usually, this the projection of the .Jacobian along the 
unstable fibers, which has direct connection with the regular part of the Perron-Frobenius operator Moreover the 
SRB measure maximizes some version of a free energy (topological pressure) : it has therefore the characteristics of 
an equilibrium state. 




(150) 



78 



References 

[1] Abeles, M. (1994). "Firing rates and well-timed events." In Domany, E., Schulten, K., and van Hemmen, J. L., editors, 

Models of Neural Networks II, chapter 3. Springer, New York. 
[2] Amit D., Gutfreund H., Sompolinsky H., "Spin glass model of neural networks", Phys. Rev. A, Vol 32, Num 2, 1007, 

(1985). Amit D., Gutfreund H., Sompolinsky H., "Storing infinite numbers of patterns in a spin glass model of neural 

networks", Phys. Rev. Letters. Vol. 55, Num. 14, 1530-1533 (1985). 
[3] Abbott, L. F., Nelson, S.B., (2000) "Synaptic plasticity: taming the beast" , Nature Neuroscience, 1178-83, Vol 3. 
[4] Abbott, L. F. and Kepler, T. B. (1990). "Model neurons: from Hodgkin-Huxley to Hopfield." In Garrido, L., editor. 

Statistical Mechanics of Neural Networks. Springer, Berlin. 
[5] Abbott, L. F. and van Vreeswijk, C. (1993). "Asynchronous states in a network of pulse-coupled oscillators". Phys. Rev. 

E, 48:1483-1490. 

[6] Adrian, E. D. (1926). "The impulses produced by sensory nerve endings". ,1. Physiol. (Lond.), 61:49-72. 
[7] Adrian, E. D. (1928). "The basis of sensation" . W. W. Norton, New York. 

[8] Amari S. (1972), "Characteristics of Random Nets of Analog Neuron-like Elements", IEEE Trans. Syst. Man. Cyb., 
SMC-2, 5. 

[9] Amari S., (1974), "A method of statistical neurodynamics", Kybernetik. 
[10] Amari S., Yoshida K., Kanatani K., (1977), "A mathematical fundation for statistical neurodynamics", SIAM Journal of 
appl. Math. 

[11] Amit D.J. "Modelling brain functions: the world of attractor neural networks", Cambridge University Press (1989). 
[12] Arnold V., "Equations diffcrcntielles ordinaires". Editions Mir, Moscou. 

[13] Arnold V., "Chapitre supplement aire de la theorie des equations differentielles ordinaires". Editions Mir, Moscou. 
[14] Arnold V., Avez A., (1967) "Problemes ergodiques de la mecanique classique", Gauthier-Vilars. 

[15] Babloyantz A., Nicolis C, Salazar J.M. [1985] "Evidence of chaotic dynamics of brain activity during the sleep cycle", 
Phys. Lett. lllA, 152-156. 

[16] Babloyantz A., Destexhe A., "Low-dimensional chaos in instance of epilepsy", Proc. Nat. Acad. Sci., USA, (1986), 3513- 
3517. 

[17] Babloyantz A., Destexhe A., "Chaos in Neural Networks", Proc. IEEE, first international conference on Neural Networks, 

M. Candill and C. Butler (eds), 4 (1987), 31-40. 
[18] Babloyantz A., Destexhe A., "The Crentzfel- Jacob disease in the hierarchy of chaotic attractors", in "From chemical to 

biological organization", M. Markus, S. MuUer, G. Nicolis. (eds). Springer series in Synergetics, 39, (1988), 307-316. 
[19] Bak P. "How Nature works: the science of Scld-Organized Criticality." , Springer- Verlag (1996), Oxford University Press 
(1997). 

[20] Chialvo D.H., Bak P., Neuroscience, 90, 1137, (1999). 

[21] Berry H., Quoy M., "Structure and dynamics of random recurrent neural networks", submitted (2005). 
[22] Bi GQ, Poo MM., J Neurosci. 1998, 18(24): 10464-72. 

[23] Basti G., Perrone A. [1989] "On the cognitive function of deterministic chaos in neural networks", IEEE, I, 657-663. 
[24] Bcnaim M., (1992), "Dynamiques d'activation et dynamiques d'apprentissage des reseaux de neurones". These de doctorat, 

Toulouse. 

[25] Bialek, W., Rieke, F., de Ruyter van Stevenick, R. R., and Warland, D. (1991). "Reading a neural code". Science, 
252:1854-1857. 

[26] Binder K., Young A., "Spin glasses : experimental facts, theoretical concepts and open questions.". Review of Modern 

Physics, Vol 58, N° 4, (1986) 
[27] Boltzmann L., "Lectures on Gas theory", Dover Publishing Co, New York, (1995). Translation by S. Brush. 
[28] Bray A. J., Moore M.A., J. Phys. C, 13 (1980) L469-476. 

[29] A Roxin, N Brunei and D Hansel (2005), "The role of delays in shaping the spatio-temporal dynamics of neuronal activity 

in large networks" , Physical Review Letters, in press 
[30] Carr (1981), "Applications of center manifold theory". Springer- Verlag: New- York, Heidelberg, Berlin. 
[31] see e.g. the web site http://elcgans.swmed.edu/ and references therein. 

[32] Cessac B., Doyon B., Quoy M., Samuelides M., (1994), "Mean-field equations, bifurcation map and route to chaos in 

discrete time neural networks", Physica D, 74, 24-44- 
[33] Cessac B., (1994), "Absolute Stability criteria for random asymmetric neural networks", J. of Physics A, 27, L927-L930. 
[34] Cessac B.,(1994), "Occurence of chaos and AT line in random neural networks", Europhys. Let, 26 (8), 577-582. 
[35] Cessac B.,"Proprietes statistiques des dynamiques de reseaux neuromimetiques". These Universite Paul Sabatier, Toulouse 

(1994). 

[36] Cessac B., "Increase in complexity in random neural networks", J. de Physique I (France), 5, 409-432, (1995). 
[37] Cessac B., Sepulchre J. A., (2004), "Stable resonances and signal propagation in a chaotic network of coupled units", 
Phys. Rev. E 70, 056111. 

[38] Cessac B., Sepulchre J. A., (2006) "Transmitting a signal by amplitude modulation in a chaotic network" Chaos 16, 
013104. 

[39] Cessac B., Sepulchre J. A., "Linear Response in a class of simple systems far from equilibrium", submitted to Physica D. 
[40] Cessac B., "Dynamical and topological effects of hebbian learning in a simple neural network model.", in preparation. 
[41] Carpenter (1977). 



79 



[42] Changcux J. P., Dchacnc S. [1989] "Neuronal models of cognitive function", Cognition 33, 63-109. 

[43] Chillingworth D.R.J. (1976), "Differentiable toplogy with a view to applications". Pitman, London. 

[44] Cohen M.A., Grossberg S., (1983), "Absolute stability of global pattern formation and parallel formation and parallel 

memory storage by competitive Neural Networks", IEEE Transaction on system, man and cybernetics, Vol SMC-13. 
[45] Cosnard M., Demongeot J., Lausberg K., Lott K., "Attractors, confincrs, and fractal dimensions: applications in neuro- 

modelling", in "mathematics applied to biology and medicine", Wucrz publishing ltd, (1993). 
[46] Crisanti A., Sompolinsky H., "Dynamics of spin systems with randomly asymmetric bonds: Langevin dynamics and 

spherical model.", Phys. Rev. A, Vol. 36, Num. 10, (1987), 4922. 
[47] Crisanti A., Sommers H.J., Sompolinsky H., "Chaos in Neural Networks: Chaotic solutions", preprint, (1990). 
[48] Cronin, J. (1987). "Mathematical aspects of Hodgkin-Huxley theory". Cambridge University Press, Cambridge. 
[49] David O., Priston K.J., "A neural mass model for MEG/EEG coupling and neuronal dynamics" (2003), Neurolmage 20, 

1743-1755. 

[50] Dauce E., Quoy M., Cessac B., Doyon B. and Samuelides M., (1998), "Self- Organization and Dynamics reduction in 
recurrent networks: stimulus presentation and learning", Neural Networks, (11), 521-533. 

[51] E. Dauce, "Adaptation dynamique et apprentissage dans des reseaux de neurones recurrents aleatoires", these troisieme 
cycle, Toulouse, (2000). 

[52] A. Guillot et E. Dauce (eds), (2002) "Approche dynamique de la cognition artificielle" , Lavoisier, Paris.. 

[53] Do Almeida A.M.O., Thouless D.J., "Stability of the Sherrington-Kirkpatrick solution of a spin glass model", J. Phys. 

A., Vol 11, Num 5, (1978), 983-990. 
[54] Doyon B., Cessac B., Quoy M., Samuelides M.,(1994), "On bifurcations and chaos in random neural networks", Acta 

Biotheoretica., Vol. 42, Num. 2/3, 215-225. 
[55] Doi S., Kumagai S., (2001),"Non linear dynamics of small scale biophysical neural networks", in Biophysical Neural 

Networks, R.R. Poznanski, ed., Mary Ann Liebert, Inc., Larchmont, NY, 261-301. 
[56] Y. Dudai, (1989), "The Neurobiology of Memory". Oxford University Press, Oxford. 

[57] Doyon B., Cessac B., Quoy M., Samuelides M., (1993), "Chaos in Neural Networks With Random Connectivity", Inter- 
national Journal Of Bifurcation and Chaos, Vol. 3, Num. 2, 279-291. 

[58] Eckmann J. P., Rucllc D., "Ergodic Theory of Strange attractors" Rev. of Mod. Physics, 57,(1985), 617. 

[59] Eckhorn R., Bauer R., Jordan W., Brosch M., Kruse W., Munk M., Reitboeck H.J. [1988] "Coherent oscillations: A 
mechanism of feature linking in the visual cortex? Multiple electrode and correlation analysis in the cat" , Biol. Cybernet. 
60, 121-130. 

[60] Edclman A., "The circular law and the probability that a random matrix has k real eigenvalues" (1993) 1-20. 
[61] Edwards, Anderson (1978), "Theory of spin-glasses" , J. Phys. A, 11, 983. 

[62] Ermentrout, G. B. and Kopell, N. (1984). "Frequency plateaus in a chain of weakly coupled oscillators". SIAM J. Math. 
Anal., 15:215-237. 

[63] Freeman W.J. "Simulation of Chaotic EEG Pattern with a dynamic model of the olfactory system", BioZ. Cyber., 56, 
(1987) 139-150. 

[64] Freeman W.J., Yao Y., Burke B.,. "Central pattern generating and recognizing in olfactory bulb: a correlation learning 

rule". Neural Networks, 1, (1988), 277-288. 
[65] FitzHugh, R. (1961). "Impulses and physiological states in models of nerve membrane" .Biophys. J., 1:445-466. 
[66] Fitzsimonds RM, Song HJ, Poo MM., "Propagation of activity dependent synaptic depression in simple neural networks" , 

Nature. 1997 Jul 31;388(6641):427-8. 
[67] Bovicr A., Gayrard V., "Rigorous results on the thermodynamic limit of the diluted Hopfield model", J. Stat. Phys., 69: 

597-627. (1993). 

[68] Gallcz D., Babloyantz A. [1991] "Predictability of human EEG: a dynamical approach", Biol. Cybern. 64, 381-392. 

[69] Gambaudo J.M., Tresser C, (1988), "Transition vers le chaos pour des applications continues de degre un du cercle", in 

"Le chaos, theorie et experiences". Collection CEA. 
[70] Geman S., "A limit theorem for the norm of random matrices", Ann. Prob., 8, (1980), 252-261. 
[71] The Genesis simulator. http://www.genesis-sim.org/GENESIS/genesis.html. 

[72] Gerstnor W.,, Kistler W.M., (2002) "Spiking Neuron Models. Single Neurons, Populations, Plasticity" Cambridge Uni- 
versity Press. 

[73] Girko V. L., "Circular law", Theor. Prob. Appl, 29, (1984), 694-706. 

[74] Gouze, J.-L. (1998). "Positive and negative circuits in dynamical systems". Journal Biol. Syst., 6(1):11-15. 

[75] Gray CM., Singer W. (1987) "Stimulus-dependent neuronal oscillations in the cat visual cortex area 17", 2nd IBRO- 

Congress, Neurosci. Suppl, 1301P. 

[76] Gray CM., Kocnig P., Engol A.K., Singer W. (1989) "Oscillatory responses in cat visual cortex exhibit intercolumnar 
synchronisation which reflects global stimulus properties " , Nature 338, 334-337. 

[77] Grimbert F. Faugeras O. (2006a) "Analysis of Jansen's model of a single cortical column", to appear in Neural Compu- 
tation. 

[78] Guckenheimer J., Holmes P., (1983), "Non linear oscillations, dynamical systems, and bifurcation of vector fields.". 
Springer- Verlag. 

[79] J. Guckenheimer, I. S. Labouriau, (1993), "Bifurcation of the Hodgkin-Huxley equations: A new twist", Bull. Math. Biol., 

55 , pp. 937-952. 
[80] J. Guckenheimer, Oliva, . 



80 



[81] J. Guckciihcinicr, Worfolk P., (1993), "Dyiianiical systems: Some computational problems.", NATO ASl, Bifurcations 
and Periodic Orbits of Vector Fields, conference proceedings and "http://arxiv.org/abs/chao-dyn/9304010". 

[82] Hassard B., (1978), "Bifurcation of periodic solutions of the Hodgkin-Huxley model for the squid geant axon", J. Theoret. 
Biol., 71, 401-420. 

[83] Hastings. (1975) 

[84] D. O. Hebb, (1949), "The Organization of Behaviour". John Wiley & Sons, New York. 

[85] Hille, B. (1992). "Ionic channels of excitable membranes". Sinauer Associates, Sunderland, Mass., 2nd edition. 

[86] Lyle J. Graham, "The Surf-Hippo Neuron Simulation System", http://www.neurophys.biomedicale.univ-paris5.fr/"' 

graham/surf-hippo- files/Surf- Hippo. README. html 
[87] Hirsch M. W. (1989) "Convergent activation dynamics in continuous time networks". Neural Networks 2, 331-349. 
[88] Hodgkin, A. L. and Huxley, A. F. (1952). "Current carried by sodium and potassium ions through the membrane of the 

giant axon of Loligo.". J. Physiol. (Lond.), 116: 449-472. 
[89] Hodgkin, A. L. and Huxley, A. F. (1952). "A quantitative description of ion currents and its applications to conduction 

and excitation in nerve membranes.". J. Physiol. (Lond.), 117:500-544. 
[90] Hopfield J. J., "Neural Networks and physical systems with emergent collective computational abilities", Proc. Natl. Acad. 

Set., USA, Vol 79, (1981), 2554-2558. 
[91] Hopfield, J. J. (1995). "Pattern recognition computation using action potential timing for stimulus representation". 

Nature, 376:33-36. 

[92] Hopfield, J. ,]., Tank, "Neural computation of decisions in optimization problems", Biol. Cybern., 52, 141-152, (1985). 
[93] Hoppensteadt F.C. and Izhikevich E.M. (1997) Weakly Connected Neural Networks, Springer- Verlag, New York 
[94] looss G. Chenciner A., (1979), "Bifurcations de tores invariants", ARMA, 69, 109-198. 
[95] Illiashenko - Weigu Li 

[96] Izhikevich E.M., "Bifurcations in brain dynamics", (1996), Ph.D., Thesis, Department of Mathematics, Michigan State 
University. 

[97] Jack J. J., Noble D., Tsien R.W., (1975) "Electric current flow in excitable cells", Oxford: Clarendon press. 

[98] Jansen B.H. Rit G. (1995), "Electroencephalogramm and visual evoked potential generation in a mathematical model of 

coupled cortical columns", Biol. Cybern., 73: 357-366. 
[99] Kolmogorov. 

[100] Katok A., Hasselblatt B., (1998), "Introduction to the modern theory of dynamical systems", Kluwer. 

[101] Keener, J. and Sneyd, J. (1998) . "Mathematical Physiology" , volume 8 of Interdisciplinary Applied Mathematics. Springer, 

New York. 

[102] Kepler, T. B., Abbott, L. F., and Marder, E. (1992). "Reduction of conductance-based neuron models". Biol. Cybern., 
66:381-387. 

[103] Kelso S.R., Ganong A.H., Brown T.H. [1986] "Hebbian synapses in hippocampus", Proc. Nat. Acad. Sci. USA, 83, 
5326-5330. 

[104] Koch, C. (1999). "Biophysics of Computation". Oxford University Press, New York. 

[105] Koch, C, O. Bernander, and Douglas, R. J. (1995). "Do neurons have a voltage or a current threshold for action potential 

initiation?", J. Comput. Neurosci., 2:63-82. 
[106] Lasalle J. P. (1968), "Stability theory for ordinary differentiel equations", J. Differential equations, 4, 57-65. 
[107] Labouriau I.S., (1989), "Degenerate Hopf bifurcation and nerve impulse H", SIAM J., Math. Anal., 20, 1-12. 
[108] Lopes da Silva F.H., van Rotterdam A., Baxts P, van Heusden E. Burr W., (1976), "Model of neuronal populations, the 

basic mechanism of rhythmicity" , M.A. Corner, D.F. Swaab (eds). Progress in brain research, Elsevier, Amsterdam, 45: 

281-308. 

[109] Mac CuUogh W.S., Pitts W., (1943) "A logical calculus of the ideas immanent in nervous activity". Bulletin of Mathe- 
matical Biophysics, 5, 11 5- 133. 

[110] MacKay R.S., Tresser C, (1986), "Transition to topological chaos for circle maps", Physica 19D, 206-237. 

[Ill] Marcus CM., Westrevelt R.M., "Dynamics of Iterated Map neural networks", Phys. Rev. A, 40, 501-504, (1989) 

[112] Marsden J.E., Mac Craken M., (1976), "The Hopf bifurcation and its applications". Springer- Verlag : New York, Heidel- 
berg, Berlin. 

[113] M.V. Mascagni, A.S. Sherman, "Numerical methods for neuronal modeling", in: Christof Koch Idan Segev (Eds.), 

Methods in Neuronal Modeling, MIT Press, Cambridge, MA, 1998. 
[114] Mezard M., Parisi G., Virasoro M.A., (1987), "Spin-glass theory and beyond". World scientific Singapore. 
[115] Milnor J., "On the concept of attractor" , Com. Math. Phys., 99, 177, (1985). 

[116] Molgedey L., Schuchardt J., Schuster H.G., "Supressing chaos in neural networks by noise", Phys. Rev. Let., Vol 69, Num 
26, (1992), 3717-3719. 

[117] Morris, C. and Lecar, H. (1981). "Voltage oscillations in the barnacle giant muscle fiber". Biophys. J., pages 193-213. 

[118] Mishchenko E.F., Rozov N. Kh. (1980) "Differential equations with small parameters and relaxation oscillations", Trans- 
lated from Russian by F.M.C. Goodspeed. New York: Plenum. 

[119] Nagumo, J. S., Arimoto, S., and Yoshizawa, S. (1962). "An active pulse transmission line simulating nerve axon". Proc. 
IRE, 50:2061-2070. 

[120] Nelson, M., Rinzel, J. (1995). "The Hodgkin-Huxley model." In Bower, J. M. and Beeman, editors. The book of Genesis, 

chapter 4, pages 27-51. Springer, New York. 
[121] Parisi G., "Asymmetric neural networks and the process of learning", J. Phys. A., 19, (1988) L675-680. 
[122] Poincare H., "Oeuvres completes", Jacques Gabay. 



81 



[123] Oram, M. W., Wiener, M. C, Lcstienne, R., and Richmond, B. J. (1999). "Stochastic nature of precisely timed spike 

patterns in visual system neuronal responses", J. NeurophysioL, 81:3021-3033. 
[124] Pollicott M., Invent. Math., 81, 413-426 (1985), Ruelle D., J. Differential Geometry, 25 (1987), 99-116. 
[125] Rieke, F., Warland, D., de Ruyter van Steveninck, R., and Bialek, W. (1996). "Spikes - Exploring the neural code". MIT 

Press, Cambridge, MA. 

[126] Rinzel, .]. (1985). "Excitation dynamics: insights from simplified membrane models" .Federation Proc, 44:2944-2946. 
[127] Rinzel, J. and Ermentrout, G. B. (1989). "Analysis of neuronal excitability and oscillations." In Koch, C. and Segev, I., 

editors. Methods in neuronal modeling. MIT Press, Cambridge, MA 
[128] Rinzel J., Miller R., (1980), "Numerical calculations of stable and unstable periodic solutions to the Hodgkin-Huxley 

equations", Math. Biosci., 49, 22-59. 
[129] Ruelle D., (1989), "Elements of differentiable dynamics and bifurcation theory". Academic Press. 
[130] Ruelle D., Takens F., "On the nature of turbulence" , Comm. Math. Phys., 20, (1971), 167-192. 

[131] Ruelle D. "Smooth dynamics and new theoretical ideas in nonequilibrium statistical mechanics" Jour. Stat. Phys. 95, 
(1999), 393-468. 

[132] Samuelides M., Doyon B., Cessac B., Quoy M.,(1996), "Spontaneous dynamics and associative learning in an asymmetric 

recurrent neural network", Math, of Neural Networks, 312-317. 
[133] Olivier Moynot, Manuel Samuelides: Large deviations and mean-field theory for asymmetric random recurrent neural 

networks Probab Theory Relat Fields 123 (2002) 1, 41-75 
[134] Quoy M., " Apprentissage dans les reseaux neuromimetiques a dynamique de base chaotique", These (1994), ENSAE, 

Toulouse. 

[135] Quoy M., Doyon B., Samuelides M., "Hebbian learning in discrete time chaotic neural networks", WCNN, Washington 
DC, (1995). 

[136] Quoy M., E. Dauce E., "Visual and motor learning using a chaotic recurrent neural network: application to the control 

of a mobile robot" , Neural Computation 2000, Berlin. 
[137] E. Dauce, M. Quoy, "Random recurrent neural networks for autonomous system design", SAB 2000, Paris, France 
[138] E. Dauce, M. Quoy, B. Doyon," Resonant spatio-temporal learning in large random neural networks", Biological Cyber- 

netics,87: 185-198,2002 

[139] Skarda C.A., Freeman W.J. [1987] "How brains make chaos in order to make sense of the world", Behav. Brain Sci. 10, 
161-195. 

[140] Shadlen, M. N. and Newsome, W. T. (1994), "Noise, neural codes and cortical organization". Curr. Opin. NeurobioL, 

4:569-579. 

[141] Sherrington D., "An introduction and overview is given of the theory of spin glasses and its application", cond- 
mat/9806289. 

[142] Soule C, "Graphic requirements for multistationnarity" , ComPlexUs, (2003);1; 123-133. 

[143] Sinai Ya. G. , "Gibbs measures in ergodic theory", Russ. Math. Surveys, 27 No 4, 21-69, (1972); Ruelle D. "Thermo- 
dynamic formalism" (1978). Reading, Massachusetts: Addison- Wesley; Bowen R. "Equilibrium states and the ergodic 
theory of Anosov diffeomorphisms" , Lect. Notes. in Math.,, Vol. 470, Berlin: Springer- Verlag (1975). 

[144] Sherrington D., Kirkpatrick S., "Solvable model of spin glass", Phys. Rev. Let. Vol 35, Nurn 26 (1975), 1792. 

[145] Smale S. (1976), "On the differential equations of species in competition". Journal of Mathematics and Biology, 3, 5-7. 

[146] Smith H.L., "Systems of ordinary differential equations which generate order preserving flows: a survey of results", SIAM 
Review, 30, 1, 87-113. 

[147] Softky, W. R. (1995). "Simple codes versus efficient codes". Curr. Opin. NeurobioL, 5:239-247. 

[148] Sompolinsky H., Crisanti A., Sommers H.J. "Chaos in random neuralnetworks" , Phys. Rev. Lett. 61, 259-262,{19S8). 
[149] Thomas R., "On the relation between the logical structure of systems and their ability to generate multiple steady 

states or sustained oscillations" , in "Numerical methods in the study of critical phenomena" , Vol 9 of Springer- Verlag in 

Synergetics, 180-193, (1981). 
[150] Thouless D.J., Anderson P.W., Palmer R. J., Phylos. Mag., 35, (1977), 593-601. 

[151] Thorpe, S., Fize, D., and Marlot, C. (1996). "Speed of processing in the human visual system". Nature, 381:520-522. 
[152] Viana M., "Stochastic dynamics of deterministic systems". 

[153] Boffeta G., Lacorata G., Musacchio S., Vulpiani A., Chaos, vol.13, pag.806, (2003) and references therein. 
[154] Williams R.F., "Expanding attractors", Publ. Math. IHES, 43; 169, (1974). 

[155] Yoshioka M., "Chaos synchronization in gap-junction-coupled neurons", ArXiv nlin. CD/0505054 (2005) 



